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ABSTRACT 


Travel  time  of  an  acoustic  signal  from  transmitter  to  receiver  provides  a 
great  deal  of  information  about  the  ocean  environment.  Variations  in  the  travel 
time  of  the  signal  may  be  caused  by  the  changes  in  the  sound  speed  along  the 
path.  Since  sound  speed  is  a  function  of  pressure,  temperature  and  salinity, 
measurement  of  this  parameter  in  acoustic  tomography  provides  a  means  to 
observe  ocean  fluctuations  through  the  use  of  inverse  techniques.  The  upcoming 
Heard  Island  Experiment  will  attempt  to  determine  the  feasibility  of  measuring 
global  warming  by  measuring  changes  in  signal  travel  time  that  may  be  caused 
by  temperature  changes  in  the  world's  oceans.  The  signals  to  be  transmitted  in 
this  experiment  are  phase  encoded  maximal-length  sequences  of  various  lengths 
which  are  well  suited  to  measurement  of  travel  time.  The  objectives  of  this 
thesis  are  to  provide  a  software  package,  in  C,  that  will  allow  participation  as  a 
receiver  in  this  experiment,  and  to  provide  a  general  capability  to  process  any 
maximal-length  sequence,  transmitted  at  any  carrier  frequency  and  with  any 
reasonable  Doppler.  A  background  on  wave  propagation,  maximal-length 
sequences,  and  Doppler  processing  are  presented  in  this  thesis. 
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I.  INTRODUCTION 


A.  THESIS  SUMMARY 

The  objective  of  this  thesis  is  to  develop  a  program  in  the  C  programming 
language  [Ref.  1]  that  can  process  acoustic  signals  used  in  ocean  acoustic 
tomography.  The  goal  of  tomography  signal  processing  is  the  precise 
measurement  of  acoustic  travel  time,  the  integral  along  the  raypath  of  inverse 
sound  speed.  Travel  time  is  a  fairly  well  understood  function  of  temperature, 
salinity,  and  pressure.  Once  variations  in  the  travel  time  are  known,  as  well  as 
the  travel  times  for  multiple  arrival  paths,  the  fluctuations  of  the  ocean  can  be 
determined  from  mathematical  inverse  techniques  [Ref.  2]. 

One  method  in  which  travel  time  can  be  measured  is  through  the  use  of 
explosive  and  implosive  devices.  These  crude  tools,  however,  are  not  exactly 
repeatable.  A  better  method  that  has  been  found  employs  the  use  of  maximal  - 
length  sequences  as  a  phase  modulation  for  bandpass  signals.  Maximal-length 
sequences  (m-sequences)  or  pseudo-random  noise  are  well  suited  for  this 
application  because  of  their  deterministic  nature,  correlation  properties,  and 
simplicity. 

The  goal  of  the  work  described  in  this  thesis  was  to  be  able  to  measure 
arrival  time  by  performing  replica  cross-correlation  of  a  transmitted  acoustic 
signal  that  has  been  p!  ase-encoded  using  m-sequences,  and  to  be  able  to  detect 
that  signal  over  any  reasonable  Doppler  shift.  Specifically,  the  programming 
package  should  be  able  to: 
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1.  Detect  any  phase  encoded  m-sequence  transmitted  with  any  carrier 
frequency. 

2 .  Scan  over  any  reasonable  Doppler  range  at  any  Doppler  bin 
spacing. 

3 .  Provide  filtering  capability  with  any  type  filter. 

4.  Employ  fast  algorithms  to  minimize  processing  time. 

5 .  Be  well  documented  and  portable  for  future  transfer  to  a  dedicated 
machine  for  real  time  processing. 

The  body  of  this  thesis  is  structured  as  follows: 

Chapter  II.  Acoustic  Wave  Propagation. 

Chapter  III.  M-Sequences  and  Hadamard  Processing. 

Chapter  IV.  Time-Varying  Doppler  Processing. 

Chapter  V.  Results  and  Conclusions. 

Chapter  II  presents  an  introduction  to  acoustic  wave  propagation  in  the  ocean 
environment.  It  includes  plane  wave  propagation,  spherical  wave  propagation, 
the  dependencies  of  sound  speed  on  temperature  and  other  parameters.  Acoustic 
raypath  determination  is  discussed. 

The  third  chapter  is  an  introduction  to  general  shift  registers  and  m- 
sequences.  The  Fast  Hadamard  Transform  is  reviewed  to  provide  the 
background  necessary  to  understand  the  programming  package. 

Chapter  IV  presents  the  signal  processing  aspects  of  the  thesis.  It  first 
presents  basic  phase  encoding  using  m-sequences,  and  then  develops  the  treatment 
of  a  signal  with  zero  Doppler.  And  finally,  the  nonzero  Doppler  case  is 
presented. 

Appendix  A  contains  the  results  from  all  Doppler  considerations.  Appendix 
B  is  the  source  code  for  the  thesis  work  as  well  as  some  supplementary 


programs.  Supplementary  programs  include:  plotting  routines  using  VAX  NCAR 
Graphics  [Ref.  3],  a  routine  for  generating  shift  register  states  and  corresponding 
m-sequences  in  a  (1,-1)  format  for  viewing  before  processing,  if  desired. 
Finally,  a  MATLAB  program  that  can  be  used  for  generating  filter  coefficients 
for  a  Butterworth  bandpass  filter  and  a  Chebychev  lowpass  filter.  The  filter 
coefficients  are  stored  to  a  file  that  can  be  read  by  the  main  program  SEQREM. 

Documentation  of  the  source  code  should  be  sufficient  for  the  typical  user  to 
understand  the  general  operation  of  the  program  without  much  difficulty,  and  is 
built  using  standard  G  functions  and  coding  to  minimize  any  portability 
problems. 

B.  THE  HEARD  ISLAND  EXPERIMENT 

The  first  direct  application  of  this  thesis  is  planned  for  the  beginning  of 
1991.  From  January  23  -  February  4,  1991  an  experiment  to  determine  the 
feasibility  of  measuring  global  warming  using  underwater  acoustic  signals  is  to 
take  place  and  is  currently  designated  the  Heard  Island  Experiment  [Ref.  4].  It  is 
hoped  that  changes  in  the  temperature  of  the  Earth's  ecosystem  can  be  measured 
by  observing  changes  in  travel  time  of  an  acoustic  signal  over  distances  which 
include  one  or  more  oceans.  Travel  time  is  related  to  water  temperature,  as  well 
as  other  parameters,  so  that  a  change  in  travel  time  might  be  related  to  warming 
of  the  Earth's  atmosphere,  which  would  also  cause  a  overall  warming  of  the 
world's  oceans.  If  these  two  parameters  can  be  related,  then  existence  of  global 
warming  might  be  determined. 

The  signals  will  be  transmitted  from  the  vicinity  of  Heard  Island  in  the 
southern  Indian  Ocean  with  raypaths  that  reach  several  receiving  sites  in  the 
Indian  Ocean,  the  Atlantic,  and  the  Pacific  as  shown  in  Figure  1.1.  The 


transmitting  vessel  will  be  the  m/v  Cory  Chouest.  Preliminary  ray  traces  have 
shown  that  one  possible  reception  site  is  the  Monterey  Peninsula  [Ref.  5]. 
Participation  is  intended  by  the  Naval  Postgraduate  School,  and  is  one  motivation 
for  this  thesis.  Of  the  many  signals  that  will  be  transmitted,  one  is  an  m- 
sequence  of  255  digits  with  a  Q  of  5,  and  a  carrier  frequency  of  57  Hz,  where  Q 
is  the  ratio  of  digit  frequency  to  carrier  frequency  and  is  normally  selected  to  be 
an  integer  number.  A  carrier  of  57  Hz  was  chosen  because  of  its  long 
propagation  distance  and  also  to  help  distinguish  it  from  the  common  frequencies 
50  and  60  Hz  that  are  generated  by  a  large  number  of  power  plants  and 
machinery  world  wide. 

Detailed  signal  specifications  are  contained  in  Reference  4,  but  some  of  the 
more  general  data  are: 

1.  HLF4LL  source. 

2 .  Output  power:  209  dB. 

3 .  Signal  to  be  transmitted  using  m-sequences 
s(t)  =  Acos(27cfg  +  M(t)y) 

where  fc  =  57  Hz. 

M(t)  is  the  m-sequence. 
y  =  45  degrees,  is  the  modulation  angle. 

4 .  Maximum  of  3-5  knots  to  maintain  way. 


4 


0  0  0  0  0  9  0 

o  o  p  <=>  o  o  P 

o  "y  Cv  c\j  T 


Figure  1.1:  Heard  Island  Raypaths. 
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II.  ACOUSTIC  WAVE  PROPAGATION 
A.  THE  WAVE  EQUATIONS 

Ocean  acoustic  signals  frequently  travel  from  source  to  receiver  along 
multiple  paths.  These  multiple  "raypaths"  are  formed  as  a  direct  result  of  the 
nonhomogeneous  structure  of  the  ocean.  However,  understanding  how  waves 
propagate  in  the  ocean  can  be  facilitated  by  first  treating  it  as  a  homogeneous 
medium  and  developing  equations  for  wave  propagation,  and  next,  by  the 
introduction  of  Snell's  Law  and  some  basic  rules  for  determining  the  raypaths  in 
a  nonhomogeneous  medium. 

Sound  travels  as  a  result  of  the  displacement  of  particles  at  some  source, 
which  causes  a  local  change  in  pressure.  The  resulting  pressure  change  then 
propagates  away  from  the  source  in  a  spherical  manner  as  a  pressure  wave  with 
velocity  c.  At  distances  sufficiently  far  from  the  source  the  wave  can  be  treated 
as  planar  over  a  finite  region. 

1.  Plane  Wave  Equation 

Assuming  the  propagation  direction  to  be  in  the  x  direction  with  the 
wave  front  formed  in  the  y-z  plane,  the  pressure  change  and  particle  speed  can 
be  related  by  the  equation  [Refs.  6,7,8] 

dp(x,t)  _  9u(x,t) 

3x  P  3t  ’ 

where  p(x,t)  is  pressure,  u(x,t)  is  particle  speed,  and  p  is  the  density  of  water. 
The  rate  of  change  of  u(x,t)  with  respect  to  x  is  given  by  [Ref.  6] 
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(2.2) 


0u(x,t)  _  _  1  8p(x,t) 

8x  B  3t  ’ 

where  B  is  the  bulk  modulus  of  elasticity.  Combining  Eq.  (2.1)  and  Eq.  (2.2) 


gives 


d  p(x,t)  B  d  p(x,t) 


(2.3) 


The  sound  speed  in  water,  c,  is  related  to  B  and  p  and  is  given  by 


2  B 
c  =  — , 


(2.4) 


Substituting  into  Eq.  (2.3)  gives  the  acoustics  plane  wave  equation  [Refs.  6,7] 


2  2 
d  p(x,t)  _  id  p(x,t) 

2  C  2~ 
0t  dx 


(2.5) 


2.  Spherical  Wave  Equation 

The  next  step  in  the  process  is  to  work  backward  from  Eq.  (2.5)  and 
find  an  expression  for  the  traveling  wave  in  spherical  coordinates.  First,  assume 
a  point  source  with  the  pressure  wave  expanding  radially.  Instead  of  the  single 
dimension  x,  Eq.  (2.5)  is  three  dimensional.  Changing  coordinate  systems  from 
linear  to  Cartesian  and  taking  the  Laplacian  in  place  of  the  partial  derivatives 
with  respect  to  position  Eq.  (2.5)  takes  the  form  [Refs.  6,7] 


(2.6) 


Assuming  no  losses,  a  uniform  traveling  wave  is  simply  a  function  of 
distance  r  from  the  source  and  time,  and  is  independent  of  angular  displacement 
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from  the  source.  Therefore,  it  is  logical  to  convert  to  a  spherical  coordinate 
system  and  use  the  spherical  form  of  the  Laplacian  operator  in  Eq.  (2.6).  After 
some  simplification  [Ref.  6],  Eq.  (2.6)  becomes 

d  (rp(r,t))  2  d  (rp(r,t)) 

— — =c — r~>  (2-7) 

at  ar 

where  r  is  the  radial  distance  from  the  source.  This  is  the  spherical  form  of  the 
wave  equation  and  is  a  function  of  only  the  distance  from  the  origin,  and  time.  It 
is  the  same  as  Eq.  (2.5)  with  p(x,t)  replaced  by  rp(r,t). 

Given  a  complex  sinusoidal  pressure  function,  a  solution  to  Eq.  (2.7)  is 

given  by 

<2-8) 

where  pm  is  the  peak  pressure  amplitude  at  unit  distance,  and  c  is  the  nominal 
sound  speed  [Ref.  6]. 

B.  RAYPATH  DETERMINATION 

There  are  two  fundamental  reasons  that  an  acoustic  signal  from  a  single 
source  may  have  multiple  raypaths  and  corresponding  arrival  times  at  a  single 
receiver.  First,  the  ocean  is  a  nonuniform  medium.  It  varies  in  depth 
(pressure),  salinity,  and  temperature.  There  are  numerous  currents,  eddies  and 
tidal  effects.  Second,  it  is  limited  vertically  by  a  sea-air  interface  and  a  sea- 
bottom  interface.  Because  of  the  first  condition,  sound  speed  is  also  nonuniform 
and  is  a  function  of  salinity,  temperature,  and  pressure.  Fortunately,  a  great  deal 
is  known  about  the  dependencies  of  sound  speed  in  water  and  the  behavior  of 
sound  waves  at  the  sea-air  and  sea-bottom  interfaces.  Applying  this  knowledge, 
and  Snell’s  Law,  raypaths  can  be  accurately  determined. 
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2.  Sound  Speed 

Sound  speed  formulation  as  a  function  of  temperature,  salinity,  and 
pressure  has  been  determined  numerically,  and  is  approximately  [Ref.  7] 

c  =  1449.2  +  4.6T  -  0.055T  +  0.00029T3  +  (1.34  -  0.01T)(S-35)  +  1.58X10'6?  ,  (2.9) 

where  c  is  sound  speed,  T  is  temperature,  P  is  gauge  pressure  of  a  column  of 
water,  and  S  is  salinity.  For  most  applications  salinity  variations  are  small  and 
can  be  neglected.  Pressure  is  a  linear  function  of  depth.  Temperature  gradients 
and  layers  are  found  throughout  the  world's  oceans  and  are  in  general  a  function 
of  depth.  Warmer  water  tends  to  reside  near  the  surface  with  colder 
temperatures  at  deeper  depths.  A  certain  amount  of  mixing  can  also  occur  to 
form  isothermal  and  isosaline  water.  Temperature  layers  can  also  be  formed  due 
to  currents  and  eddies.  Sound  waves  can  be  thought  of  as  tending  to  bend  toward 
points  of  lower  sound  speed,  as  will  become  apparent  in  the  next  section  on 
Snell's  Law. 

2.  Snell's  Law 

A  raypath  through  the  ocean  medium  can  be  determined  by  considering 
the  sound  speed  versus  depth  to  be  a  continuous  function,  implying  a 
continuously  stratified  medium.  Changes  in  propagation  path  can  be  computed 
by  treating  the  entire  stratified  medium  as  a  set  of  n  layers,  and  determining  the 
refraction  from  medium  1  to  medium  2,  then  2  to  3,  and  so  on  to  medium  n. 
This  is  done  through  the  use  of  Snell's  Law,  which  is  derived  in  several  texts 
[Refs.  6,7].  Stated  here,  a  sound  wave  (ray)  will  be  refracted  from  medium  to 
medium  according  to  the  relation 
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> 


(2.10) 


sin(0j)  sin(02)  sin(0n) 

ci  c2  cn 


where  0]  is  the  incident  angle,  measured  from  the  ray  to  the  vertical,  and  02, 

0n  are  refraction  angles.  The  constants  ci,  C2,  ...  cnare  the  sound  speeds  in 
medium  1,  medium  2,  and  medium  n,  respectively.  It  is  assumed  that  the  energy 
in  the  ray  is  constant  through  the  boundary. 

As  an  example,  consider  the  case  where  cj  >  C2.  In  this  case,  the  sound 
speed  in  medium  1  is  faster  than  that  of  medium  2.  Rearranging  F.q.  (2.10)  in 


the  following  manner 


(2.11) 


it  can  be  seen  that  ratio  of  C2  to  ci  is  less  than  one.  For  the  equality  to  hold,  02 
must  be  a  smaller  angle  than  that  of  the  incident  angle,  0j.  This  means  the  ray  is 
"bending"  downward  towards  the  medium  with  -he  lower  sound  speed.  Similarly 
the  ray  bends  upwards  for  the  case  when  C2  >  c..  In  other  words,  sound  waves 
bend  toward  the  region  of  slower  sound  speed. 

In  reality,  when  sound  impinges  on  a  boundary  between  two  mediums 
the  sound  waves  are  both  reflected  and  transmitted.  Application  of  Snell's  Law 
remains  the  same  for  the  transmitted  wave,  and  the  reflected  wave  has  a 
reflection  angle  that  is  equal  to  the  incident  angle.  The  transmitted  angle  remains 
a  function  of  the  sound  speeds  of  the  medium*  on  eithr  sr'cle  of  the  boundary  and 
is  given  by  Eq.  (2.10).  However,  the  amplitude  of  u»e  transmitted  wave  does  not 
equal  that  of  the  incident  wave,  because  con:,,  nation  of  energy  must  apply. 
Therefore,  the  amplitude  of  the  reflected  and  transmitted  waves  must  sum  to  that 
of  the  incident  wave  at  the  boundary.  Assuming  the  density  of  medium  1  and  2 
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to  be  the  same,  the  ratio  of  the  amplitude  of  the  reflected  wave  to  that  of  the 
incident  wave  is  given  by  [Ref.  6] 


r  c2  sin(0  j)  —  Cj  sin(02) 

1  ci  sin(6  j)  -  -  .?*m(92)’ 


(2.12) 


where  R  is  the  reflected  amplitude,  i  is  tk  us  amplitude,  0i  is  the  incident 
angle.  02  is  the  transmitted  angle,  and  ci,  a?.d  -  _  are  the  sound  speeds  in  medium 
1  and  2,  respectively.  Similarly,  the  ra:  ■>  the  transmitted  to  incident 
amplitude  is  given  by 

t  2  c2  sintO,; 

T  = - ,  (2.13) 

1  C2  sin^B  j)  *t-  c j  sin(02) 


where  T  is  the  transmitted  amplitude.  For  the  case  when  the  mediums  do  not 
have  the  same  density,  the  sound  speeds  c;  in  Eq.  (2.12)  and  Eq.  (2.13)  are 
simply  replaced  by 

zi  =  Pic->  (2.14) 


where  p, ,  and  Z,  are  density  and  impedance  of  the  i*  medium,  respectively. 

For  rays  incident  at  li  „  sea  surface  the  sound  wa”e  is  considered  to  be 
totally  reflected,  and  the  amount  of  reflection  at  the  ocean  bottom  is  based  on  the 
local  bottom  consistency.  There  are  two  special  conditions  that  can  be  applied 
with  Snell's  Law.  In  the  first  case,  when  c2  >  c2  and  the  ratio  of  cj  to  c2  becomes 
very  large,  the  transmission  angle  02  must  approach  zero.  In  this  case,  the  wave 
transmits  perpendicular  to  the  boundary.  In  the  second  case,  when  c2  >  Cj,  there 
is  a  condition  when  the  transmitted  angle  must  be  negative  to  satisfy  Eq.  (2.10). 
However,  02  cannot  be  negative.  This  is  the  condition  for  total  reflection  and  the 
reflected  angle  equals  the  incident  angle,  and  B2  is  zero. 
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3.  Travel  Time 

Once  the  sound  speed  as  a  function  of  depth  has  been  determined,  the 
ray  path  can  be  determined.  Travel  time  can  then  be  evaluated  as  an  integral 
over  the  raypath.  Travel  time,  Tj,  for  a  particular  ray  i  from  point  a  to  b  is 
[Ref.  9] 


b 


ds 

c(x,y,z)  ’ 


(2.15) 


where  sound  speed  is  assumed  to  be  independent  of  time  t.  If  th  assumption  is 
made  that 


c(x,y,z)  =  c0(x,y,z)  +  6c(x,y,z) ,  q  I  ^ 

where  co(x,y,z)  is  an  assumed  background  sound  speed  field,  and  8c(x,y,z)  is  an 
unknown  sound  speed  field,  Eq.  (2.15)  can  then  be  written  as 

b 


‘'i“^'oi+5T,=  f  b — 7— — r  -  f  5t'y'^dS. 

I  C»(X’y-Z)  !  4(x,y,z) 


(2.17) 


for  5c  «  c0.  The  travel  time  T;  has  been  broken  down  into  a  travel  time  To,  cute 
solely  to  the  assumed  sound  speed  field  co,  and  8Tj,  a  linear  function  of  the 
unknown  sound  speed  field  perturbation  5c(x,y,z).  Variation  in  the  arrival  time 
T,  is  a  perturbation  of  the  arrival  time  due  to  perturbations  in  the  sound  speed 
along  the  path  By  measuring  these  variations  and  using  linear  inverse 
mathematical  techniques  it  is  possible  tc  determine  some  desired  characteristic  of 
this  geophysical  problem  [Refs.  2,9]. 
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III.  M-SEQUENCES  AND  HADAMARD  PROCESSING 


A.  INTRODUCTION 

An  important  parameter  in  underwater  acoustic  tomography  is  the  arrival 
time  of  a  transmitted  signal.  A  very  usefm  means  for  measuring  arrival  time  is 
to  apply  an  impulsive  excitation  to  the  system,  in  this  case  the  ocean 
environment.  An  impulse  can  be  generated  by  explosive  or  implosive  sources 
cheaply  and  easily,  but  suffer  from  unevenness  in  the  frequency  spectrum,  and 
repeatability.  Another  method  Involves  the  use  of  pseudorandom  noise. 
Psuedorandom  noise  is  deterministic  ana  is  formed  by  a  sequence  of  ones  and 
zeros  known  as  a  maximal-length  sequence  or  m-sequence  [Refs.  10,11].  re¬ 
sequences  are  typically  transmitted  by  using  the  ones  and  zeros  of  the  sequence 
([1,0]  mapped  to  [-1,1])  to  phase  encode  a  carrier  signal  [.'^ef's.  J  2,13,14]. 

M-sequences  are  generated  using  a  simple  binary  shift  register  that  is  formed 
uncording  to  the  desired  primitive  polynomial.  Since  the  signal  being  transmitted 
is  known,  a  simple  autocorrelation  can  be  formed  uring  a  bank  of  matched 
filters.  The  resulting  correlation  is  impulse-like,  and  has  a  shorter  duration  than 
the  originally  transmitted  signal  and  has  a  much  larger  amplitude,  which  makes 
measurement  of  arrival  time  an  easier  task. 

The  following  sections  give  a  brief  description  of  shift  register  fundamentals, 
m-sequences,  and  fast  Hadamard  processing  of  m-sequences. 
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B.  SHIFT  REGISTER  FUNDAMENTALS 

A  general  sequence  shift  register  generator  is  shown  in  Figure  3.1.  A 
sequence  of  ones  and  zeros  can  be  output  from  any  one  of  the  registers  starting  at 
any  register  state,  except  zero.  In  what  follows  all  output  sequences  are  taken 
from  the  Oth  register.  The  sequence  of  bits  generated  is  periodic  and  will  repeat 
with  some  period  L,  based  on  the  structure  of  the  generator,  as  described  by  a 
polynomial  such  as 

g(D)  =  bnDn+bn.1Dn  1+...+b1D+b0  >  (32) 

where  D  is  the  unit  delay  and  bj-  are  the  feedback  weighting  coefficients.  The 
weighting  coefficients  take  on  values  of  one  or  zero  indicating  connection  or  no 
connection  to  the  k^1  register.  All  arithmetic  operations  for  polynomials  are 
performed  using  finite-field  arithmetic  (modulo  two)  and  are  described  in  more 
detail  in  Reference  10. 


Figure  3.1:  General  sequence  shift  register  generator  for  Eq.  (3.1). 

An  example  of  a  third  order  sequence  shift  register  generator  is  shown 
in  Figure  3.2,  and  its  corresponding  generating  polynomial  is 

g(D)  =  D3  +  D  +  1  .  (3.2) 

In  this  case  Eq.  (3.2)  happens  to  be  a  primitive  polynomial,  where  a  primitive 
polynomial  is  any  polynomial  that  will  not  repeat  its  register  state  until  after 
2n-l  delays,  where  n  is  the  number  of  delays  in  the  shift  register  generator. 
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Therefore,  a  shift  register  that  is  defined  by  a  primitive  polynomial  will  generate 
a  sequence  that  is  of  maximal-length  and  is  known  as  a  maximal-length  sequence 
or  m-sequence. 


Figure  3.2:  Shift  register  realization  of  Eq.  (3.2). 

Continuing  with  this  example,  the  sequential  output  from  the  registers  of 
figure  3.2  is  then  given  as  in  Figure  3.3,  where  the  initial  register  contents  is 
arbitrarily  set  to  a2=l,  aj=0,  a0=0  [Ref.  9j. 

The  corresponding  m-sequence  is  then  obtained  by  taking  any  one  column 
from  the  register  states  top  to  bottom  (or  bottom  to  top).  Note,  that  the 
combination  of  [000]  never  occurs.  This  is  because  an  initial  value  of  zero  will 
not  allow  any  transition  in  state  and  can  be  easily  verified  by  inspection. 


Cycle 


1 

2 

3 

4 

5 

6 

7 

8 


«2 


1  0  0 
0  1  0 
0  0  1 
1  0  1 
1  1  1 
1  1  0 
0  1  1 
1  0  0 


Figure  3.3:  Shift  register  contents  when  generating  a  third  order  m- 
sequence.  The  eighth  cycle  shows  that  the  register  begins  to  repeat. 

The  characteristics  of  the  m-sequence  is  unchanged  whether  it  is  transmitted 

in  the  forward  or  reverse  direction,  but  when  performing  the  fast  Hadamard 

Transform,  reviewed  in  the  next  section,  the  direction  in  which  the  sequence  is 

transmitted  is  significant.  Therefore,  the  top  to  bottom  sequence  will  be 

designated  the  "forward"  code  and  the  bottom  to  top  the  "reverse"  code, 

forward  m  =  0011101 

reverse  m=  1011100.  (3.3) 

M-sequences  have  several  desirable  characteristics.  Its  correlation  function 
is  triangular  in  shape,  see  Figure  3.4,  and  short  in  duration.  It  is  deterministic, 
once  the  polynomial  is  determined  all  output  is  known.  It  can  be  implemented 
easily  and  is  periodic.  However,  one  major  drawback  is  that  its  autocorrelation 
requires  N  multiplies,  and  since  the  arrival  time  is  never  known,  the  input  signal 
must  be  correlated  with  all  N  shifted  versions  of  the  original  sequence  which 
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requires  N2  multiplies.  This  can  be  overcome  by  us^  of  the  fast  Hadamard 
transform  (FHT).  discussed  in  the  next  section. 


Figure  3.4:  Autocorrelation  function  for  a  m-sequence  of  length  N. 


C.  THE  FAST  HADAMARD  TRANSFORM 

The  autocorrelation  of  the  input  data  sequence  with  all  possible  shifted 
versions  of  the  original  m-sequence  can  be  written  as 


MD  = 


0  0  1  i  1  0  1 
10  0  1110 
0  10  0  111 
10  10  0  11 
110  10  0  1 
1110  10  0 
0  1110  10 


(3.4) 


where  M  is  a  matrix  whose  rows  are  the  shifted  versions  of  the  forward  code  and 
D  is  the  input  data  vector.  This  product  requires  N2  multiplications.  To  reduce 
this  number,  a  Hadamard  transform  can  be  used.  The  third  order  Hadamard 
matrix  is  given  by 
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11111111 
1-1  1-1  1-1  1-1 
1  1-1-1  1  1-1-1 
1-1-1  1  1-1-1  1 
111  l-l-l-l-l 
1-1  1-1-1  1-1  1 
1  l-l-l-l-l  1  1 
1-1-1  1-1  1  1-1 


(3.5) 


where  H  can  be  formed  recursively  from 

=  !=[!].  Hitl 


(3.(5) 


and  by  performing  a  simple  mapping  of  (1,-1)  to  (1,0).  H  can  be  written  in  a 


form  that  is  easier  for  most  people  to  work  with  [Refs.  15,16,17], 

"0  0  0  0  0  0  0  o' 
01010101 
00110011  - 
_  01100110 
00001111 
01011010 
00111100 
011010  0  1  • 


(3.7) 


It  is  can  also  be- shown  .-[Refs.  16,17]  that  the  matrix  H  can  be  factored  into 
the  product  of  two  matrices  consisting  of  a  binary  count,  in  this  case  from  0  to  7, 


0  0  0 

0  0  1 

0  10 

"o  0  0  0  1  1  1  1" 

Oil 

00110011 

10  0 

01010101 

10  1 

1  1  0  j 
111 

(3.8) 


The  matrix  M,  Eq.  (3.4),  can  also  be  factored  into  two  matrices,  L  and  $,  given 
by  [Refs.  16,17] 
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LS 


fl  0  0 
>010 
1  0  1 
-  1  0 
Oil 
111 
10  1 


0  0  1110  1 
10  0  1110 
0  10  0  111 


(3.9) 


It  is  easily  verified  that 


M  =  LS. 


(3.10) 


Note  that  by  adding  a  leading  column  of  zeros  to  S,  forming  S’,  and  a  leading 
row  of  zeros  to  L,  forming  L1.  All  possible  combinations  of  ones  and  zeros  are 
formed  in  L',  S',  as  in  A,  with  the  differences  between  the  matrices  being  the 
order  in  which  they  occur.  It  is  straight  forward  to  find  a  matrix  P  such  that 
S’=AT  P  and  another  matrix  U  such  that  L'=UA.  Combining  these  results  maps 
M'  to  the  Hadamard  matrix  as  [Ref.  16,17] 


M'  =  L'S'  =  UAAtP  =  UHP,  (3.11) 

where  M'  is  the  M  matrix  with  an  appended  column  and  row  of  zeros  in  the  first 
row  and  column.  Recall  that  the  correlation  was  performed  by  forming  the 
product  of  the  M  matrix  and  the  data  vector  D.  By  replacing  M  with  M'  and 
forming  a  new  data  vector 


D'  = 


a 

b 

c 

d 

e 

f 

g 

h 


(3.12) 


and  after  substituting  Eq.  (3.1 1)  for  M'  the  correlation  becomes 

R’  =  M'D’  =  UHPD’.  (3.13) 
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From  this  it  can  be  seen  that  the  matrix  M  is  not  needed  and  the  Hadamard 
matrix  can  be  used  in  its  place.  By  using  the  Hadamard  matrix  in  the  form  of 
Eq.  (3.5)  and  forming  a  product  with  D',  the  result  becomes  a  sum  of  the 
individual  terms  of  D’  with  +  and  -  weights  and  has  the  form 


HD’ 


a+b+c+d+e+f +g+h 
a-b+c-d+e-f +g-h 
a+b-c-d+e+f-g-h 
a-b-c+d+e-g-g+h 
a+b+c+d-e-f-g-h 
a-b+c-d-e+f-g+h 
a+b-c-d-e-f +g+h 
a-b-c+d-e+f +g-h 


(3.14) 


which  when  written  in  the  form  of  a  flow  graph  having  the  same  form  as  the 
well  known  fast  Fourier  transform  shown  in  Figure  3.5,  where  the  complex 
multiplies  are  set  to  one  [Refs.  9,16,17].  This  is  known  as  the  fast  Hadamard 
transform  (FHT). 

All  that  is  left  is  to  determine  the  permutation  matrices  P  and  U.  By  looking 
at  how  P  and  U  must  be  formed,  it  can  be  seen  that  P  must  have  ones  at  the 
indices  [Ref.  15,16] 

ROW  01234567 
COLUMN  0  2  1  4  6  7  3  5, 

and  U  must  have  ones  at  the  indices 


ROW  0  4  2  1  6  3  7  5 

COLUMN  0  1  2  3  4  5  6  7, 

which  correspond  to  the  binary  values  of  S'  and  L’.  It  turns  out  that  the  matrices 
S'  and  L'  provide  all  of  the  necessary  information.  All  that  is  needed  to  perform 
the  desired  multiplication  is  to  permute  the  input  data  vector  according  to  S', 
form  the  FHT,  and  permute  again  according  to  L’.  L  and  S  are  generated  from 
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the  primitive  polynomial  that  defines  the  m-sequence  using  the  generators  shown 
in  Figure  3.6  for  the  forward  code  and  Figure  3.7  for  the  reverse  code  [Ref.  15]. 
The  modification  to  form  S’  and  L'  simply  involves  adding  leading  zeros  to  S 
and  L.  Since  the  correlation  is  always  over  2n-l  data  points,  a  leading  zero  is 
added  to  the  input  data  vector  D  forming  D‘  adding  no  new  data.  After 
processing,  D’  holds  the  correlation  function  with  the  zero  position  containing 
the  average  level,  which  may  be  used  to  remove  the  bias  in  the  autoconeiaiion 
function  [Ref.  11]. 

Basic  Element 


a+b+c+d+e+f+g+h 
a-b+c-d+e-f+g-h 
a+b-c-d+e+f-g-h 
a-b-c+d+e-f-g+h 
a+b+c+d-e-f-g-h 
a-b+c-d-e+f-g+h 
a+b-c-d-e-f+g+h 
a-b-c+d-e+f+g-h 

Figure  3.5:  Basic  fast  Hadamard  transform  element  for  cascading 
additions  and  the  full  diagram  for  an  eight  point  FHT. 
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IV.  TIME-VARYING  DOPPLER  PROCESSING 


A.  PHASE  ENCODING  OF  M-SEQUENCES 

A  simple  method  in  which  m-sequences  can  be  transmitted  is  to  phase  encode 
the  digits  using  a  mapping  of  (0,1)  to  (1,-1)  and  appropriately  selecting  a  phase 
angle  which  will  maximize  the  signal-to-noise  performance.  The  phase  angle 
that  optimizes  this  is 

\j/  =  tan'^VN),  (4.1) 

where  \\r  is  the  phase  modulation  angle,  and  N  is  the  length  of  the  sequence  [Ref. 
11].  However,  for  the  Heard  Island  Experiment,  \(/  =  45  degrees  is  used  to 
simplify  SNR  estimation.  The  transmitted  signal  s(t)  takes  the  form 

s(t)  =  Acos(27tfct  +M(t)\j/),  (4.2) 

where  A  is  the  amplitude,  fc  is  the  earner  frequency,  and  M(t)  takes  on  values  of 

+1  or  -1  depending  on  the  m-sequence  code.  The  minimum  length  of  time  M(t) 
remains  constant  is  the  digit  duration  d. 

It  is  normal  practice  to  chose  an  integer  number  of  cycles,  Q,  of  the  carrier 
frequency  per  digit. 

B.  SIGNAL  PROCESSING  WITH  ZERO  DOPPLER 

At  the  receiver  the  signal  is  normally  sampled  at  an  integer  multiple  of  the 
carrier  frequency  fc,  i.e.  with  fs  =  mfc,  where  m  is  chosen  to  be  an  integer 

greater  than  two,  to  ensure  that  the  Nyquist  criterion  is  satisfied,  and  fs  is  the 

sampling  frequency.  Sampling  in  this  fashion  simplifies  processing  and 
interleaving  of  data,  since  there  is  an  integer  number  of  data  points  per  digit  d. 
The  received  signal  is  r(t)  =  s(t  -  T,),  assuming  no  attenuation,  no  dispersion,  and 


24 


a  single  raypath.  The  sampled  input  r(n)  is  then  den  adulated  to  remove  the 
carrier  from  the  signal.  The  received  signal  r(n)  will  arrive  at  some  unknown 
arrival  time  Ti,  and  must  be  multiplied  by  both  the  sine  and  cosine  functions  to 
form  the  in-phase  and  quadrature  components  p(n)  and  q^n),  respectively  [Ref. 
18J.  For  p(n)  this  then  gives 

Inisi  '  n  27tf~n 

P(n)  =  Acos(— 4-  +  M(f  -  Tj)\i/)cos(— t-— ) 

*s  rs  tS 

p(n)  =  |  [cos(M(£  -  Ti)V)  +  cos(— |r~  +  M(jr  -  Tj)^)^  (4  3) 


and  similarly  for  q(n) 


A  r.  47tf-n 

q(n)  =  j  [sin(M(~  -  Tj)\p)  -  sin(— 


+  M(~Ti)ii/)]. 


(4.4) 


The  high  frequency  components  are  then  removed  by  lowpass  filtering  with  the 
cutoff  frequency  greater  than  the  digit  bandwidth.  The  resulting  waveforms  are 


p(n)=ycos(M(~Ti)V) 

q(n)=|sin(M(^-Ti)V), 


(4.5) 


which  are  constant  over  the  length  of  the  digit  of  M(n). 

Because  of  the  lower  digit  bandwidth,  decimation  in  time  may  be  performed 
without  any  loss  of  information  [Ref.  19].  This  significantly  reduces  the  data 
rate  and  storage  requirements  of  the  processed  data.  Further  processing  follows 
the  scheme  of  the  previous  chapter.  The  input  data  vector  is  permuted,  the  FHT 
is  performed,  and  permuted  again  to  the  correct  order.  Multiple  samples  per 
digit  must  be  accounted  for  by  interleaving  when  performing  the  FHT  on  a  given 
sequence  length.  Finally,  magnitude  and  phase  are  computed  in  the  normal 
fashion. 
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The  autocorrelation  function,  see  Figure  3.4,  for  m-sequences  has  a  negative 
DC  level  for  all  digit  positions  except  the  first.  A  correction  may  be  made  which 
will  adjust  this  level  to  zero.  Knowing  the  phase  modulation  angle,  the  DC  bias 
can  be  removed  by  adding  the  correction 

^cor =  fi^avg  »  (4-6) 

where  Ccor  is  the  complex  level  adjustment,  Lavg  is  the  average  DC  level  from 

the  FHT,  and  fj  is  given  by  the  relation  [Ref.  11] 

f  _  .  (N+l)  tanCy)  +  j(N  -  tan2(y)) 

1  J  N2  +  tan  2(v)  ’  (4-7) 

where  N  is  the  length  of  the  m-sequence  and  is  the  phase  modulation  angle. 
Since  all  parameters  of  Eq.  (4.7)  are  known,  fj  is  computed  once  at  system 
initialization  and  Ccor  requires  one  complex  multiplication  for  each  data  length 

processed. 

A  block  diagram  showing  the  process  is  given  in  Figure  4.1.  Note  that  in 
this  case  a  Butterworth  filter  of  order  ten  was  used  in  the  passband  and  a 
Chebychev  filter  of  order  five  was  used  for  the  lowpass  filter.  Any  filter  type 
can  logically  be  substituted. 


26 


Figure  4.1:  Block  diagram  of  signal  processing  with  zero  Doppler. 
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C .  SIGNAL  PROCESSING  WITH  DOPPLER 

Generally  in  acoustic  tomography  the  receiving  or  transmitting  ship 
must  maintain  way  and/or  buoys  are  acted  on  by  the  forces  of  currents,  tides, 
wind  etc.,  all  of  which  cause  relative  motion  between  the  receiver  and 
transmitter.  The  Doppler  shift  that  is  then  imparted  on  the  signal  must  be 
corrected  for  at  the  receiver.  For  the  relatively  slow  and  uniform  speeds 
encountered  in  acoustic  tomography  it  is  sufficient  to  compensate  for  first  order 
Doppler,  and  neglect  acceleration  effects  [Refs.  20,21]. 

When  Doppler  is  present  a  transmitted  signal  s(t)  may  be  received  as 

r(0  =  s[(l  +”)  t-  Tj],  (48) 

where  v  is  the  velocity  in  meters/second,  c  is  the  nominal  sound  speed  in  water  in 
meters/second,  and  Tj  is  the  delay  from  transmitter  to  receiver  as  in  section  B, 
where  attenuation,  dispersion,  and  a  single  raypath  are  assumed.  The  magnitude 
of  the  spectrum  of  r(t)  can  be  given  as 

|R®|« 

Sampling  at  an  integer  multiple  of  the  carrier  frequency  assuming  zero  Doppler 
and  replacing  the  time  dependency  with  the  discrete  index  n  gives 


S[ — - ] 

a+v/c> 


(4.9) 


r(n)  =  s(f--Ti). 

*s 


(4.10) 


With  Doppler,  the  received  signal  is  then  given  by 
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(4.11) 


r(n)  =  s[(- 


» a + v/o) 


)'T;1 


Adjusting  the  sampling  frequency  by  the  same  Doppler  shift  gives  a  new 
sampling  frequency  of 

fs=fs(l+%)>  (4.12) 

which  forms  a  new  received  signal 

n(l+v/c) 

-)-T;l 

(4.13) 


r'(n)  =  s[(- 


f’c 


Substituting  Eq.  (4.12)  into  Eq.  (4.13)  yields 


r’(n)  =  s(t-  -Tj)  =r(n) 


(4.14) 


At  this  point  the  results  are  identical  and  processing  of  the  data  can  proceed  as  in 
the  zero  Doppler  case  [Ref.  20]. 

Two  methods  can  be  employed  to  perform  the  adjustment  to  the 
sampling  frequency  as  discussed.  The  first  is  to  use  a  bank  of  sampling  devices 
that  sample  at  a  set  of  sampling  frequencies  that  cover  the  desired  Doppler  range 
such  that  each  device  samples  at  frequencies  This  method  requires  a 

great  deal  of  hardware  and  is  inflexible  and  expensive  to  implement.  The  second 
method  is  to  frequency  shift  the  received  signal  and  interpolate  between  samples 
at  the  new  sampling  frequency  to  predict  the  value  of  the  signal  as  if  the  desired 
sampling  frequency  had  been  used  [Refs.  20,21], 

Doppler  resolution  is  1/TA  Hz,  where  TA  is  the  analysis  interval  [Ref. 
20].  For  example,  a  N=255  length  m-sequence  with  Q=5  and  fc=57  Hz  giving  a 
Ta  of  22.368  seconds  and  a  resolution  of  0.04471  Hz,  which  from  Eq.  (4.12) 
corresponds  to  a  Doppler  shift  of  +  or  -  2.3  knots.  Figure  4.2  is  an  ambiguity 
plot  derived  from  this  example.  The  plot  shows  signal  arrival  time  for  the  given 
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m-sequence  versus  signal  Doppler  over  a  +  or  -  5  knot  range.  As  can  be  seen  in 
the  plot  signal  strength  decreases  as  the  carrier  frequency  moves  away  from  the 
zero  Doppler  case  and  is  not  detectable  at  about  2  knots  on  either  side  of  center. 

Frequency  shifting  in  the  frequency  domain  is  equivalent  to  multiplying 
the  samples  or  the  demodulates  in  the  time  domain  by  exp(j2ftnfd/fs)  [Ref.  20], 

where  fd  is  the  Doppler  shift  of  the  carrier  signal  in  Hertz.  The  Doppler  shift  fd 
is  determined  directly  from  the  expected  Doppler  and  is  given  by 


(4.15) 


For  a  sampling  frequency  that  is  four  times  the  carrier  the  multiplication 
exponential  can  then  be  related  to  the  Doppler  speed,  v  (meters/second),  and  the 
multiplication  of  the  demodulates  in  the  time  domain  is  by  exp(jrtnv/2c).  The 

Aobicuity  Plot 


Figure  4.2:  Ambiguity  plot  of  Arrival  Time  vs.  Doppler  for  a 
theoretical  sequence  with  N=255,  Q=5,  and  fc=57  Hz. 
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data  is  then  inteipolated  at  the  new  sampling  frequency  corresponding  to  the 
original  sampling  frequency  shifted  by  the  same  amount  of  Doppler  as  given  by 
Eq.  (4.12). 

The  Doppler  processing  scheme  in  this  thesis  is  performed  as  follows. 
Frequency  shifting  is  performed  by  combining  the  shift  with  the  demodulation 
process.  This  reduces  the  number  of  complex  multiplies  by  combining  the 
modulation  angle  and  the  frequency  shift  via  simple  addition,  see  Figure  4.3. 
The  resampled  signal  is  then  linearly  interpolated  [Ref.  21].  Since  the  Doppler 
shift  is  unknown,  except  within  some  range,  each  data  set  is  processed  over  the 
entire  Doppler  range,  where  the  Doppler  step  size  is  set  in  knots.  Frequency 
shifting  and  then  interpolating  results  in  the  following  expression  for  the 


resampled  signal  [Ref.  21] 


s(ta+h)  =  s(ta)ej27:fdt*  +  --  [s(tb)ej2nfdtb  -  sCt^2^'*  \ 


(4.16) 


where  ta<  ta+h<tb  and  time  replaces  the  discrete  index  n  for  clarity.  The  times  ta 
and  tb  correspond  to  sampled  data  points,  and  h  is  the  time  spacing  from  ta  to  the 
point  to  be  interpolated.  The  position  h,  varies  from  interpolation  point  to 
interpolation  point,  according  to  the  sampling  period  l/fs\  as  time  progresses. 
Other  techniques  for  interpolating  between  samples  are  available  but  are  slower 
and  more  complex  and  will  not  be  discussed  here. 

At  this  point  processing  continues  as  in  the  zero  Doppler  case.  Figure 
4.3  shows  the  process  in  block  diagram  form  with  an  additional  block  for 
performing  coherent  averaging  to  increase  signal  processing  gain,  if  required  or 
desired  [Ref.  18]. 

Figure  4.4  shows  the  Doppler  resolution  for  a  transmitted  signal  with 
-1.4  knots  of  Doppler,  using  the  parameters  given  in  the  above  example  without 
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noise,  and  was  processed  by  this  scheme.  Note  that  these  results  are  plotted  on  a 
different  scaling  than  that  of  Figure  4.2  but  show  the  same  results.  The  Doppler 
resolution  is  again  approximately  +  or  -  2  knots  as  expected. 


Figure  4.3:  Block  diagram  of  signal  processing  with  Doppler  and 
coherent  averaging  incorporated. 
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ARRIVAL  TIME  VS  DOPPLER 
-5  TO  5  KNOTS  DOPPLER 


Figure  4.4:  Ambiguity  plot  for  Arrival  Time  vs.  Doppler  for 
sequence  length  N=255,  Q=5,  and  fc=57  Hz,  processed  by  SEQREM 

program  with  no  signal  noise. 


33 


V.  RESULTS  AND  CONCLUSIONS 


A.  RESULTS 

The  objective  of  this  thesis  was  to  develop  a  program  in  the  C  programming 
language  that  can  process  any  m-sequence  at  any  carrier  frequency  over  any 
acceptable  Doppler  range.  This  objective  has  been  fully  met.  The  package 
SEQREM  has  been  developed,  which  will  scan  over  any  specified  Doppler  range 
and  specified  Doppler  step  size.  In  cases  where  no  Doppler  is  expected  the 
program  skips  the  Doppler  scanning  procedure,  reducing  run  time. 
Programming  allows  for  any  carrier  frequency,  any  type  filter  of  any  order,  up 
to  twenty,  and  any  size  m-sequence,  with  dynamic  memory  space  allocation  for 
array  manipulation,  limited  only  by  the  available  computer  memory.  It  permits 
coherent  averaging  of  sequence  segments  as  specified  at  initialization,  improving 
processing  gain,  and  computes  and  removes  the  DC  bias  in  the  autocorrelation 
function.  Decimation  in  time  is  performed  which  significantly  reduces  off-line 
storage  requirements  for  processed  data. 

Appendix  A  contains  a  complete  set  of  processing  results  for  all  three  cases, 
closing,  opening,  and  zero  Doppler.  In  all  of  the  plots  a  255  digit  sequence  is 
used  with  Q  =  5,  and  a  carrier  of  57  Hz,  corresponding  to  one  of  the  signals  to 
be  transmitted  in  the  Heard  Island  Experiment.  The  sampling  frequency  ;s  four 
times  the  carrier  frequency  or  228  Hz.  The  simulated  data  assumes  white 
Gaussian  noise.  The  plots  are  labeled  with  the  Doppler  that  was  processed  and 
with  the  sequence  period,  adjusted  for  the  expansion  or  compression,  according 
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to  the  Doppler  bin.  Decimation  in  time  was  performed  at  a  ratio  of  ten  to  one. 
Peaks  correspond  to  the  center  of  detection  of  the  signal  arrival  within  the 
period.  The  leading  edge  of  the  peak  is  the  actual  arrival  time  and  can  be 
extracted  by  an  edge  detection  program.  Sequential  frames  are  aligned  one 
behind  the  other  and  represent  increasing  time  along  the  y  direction.  All  plots 
are  in  magnitude.  Phase  plots  are  omitted  for  brevity.  Some  specific  results  are 
contained  here. 

Figure  5.1  shows  the  detection  of  a  zero  Doppler  signal  that  has  been  shifted 
in  time  and  has  some  arbitrary  arrival  phase.  The  signal  is  shifted  approximately 
three  digits  with  an  SNR  =  0  dB. 

The  sequence  removal  process  is  linear  and  can  detect  multiple  signal 
arrivals.  Figure  5.2  shows  this  clearly  for  two  arrivals  spaced  27  digits  apart 
(approximately  6  seconds)  with  SNR  =  -15  dB.  Figure  5.3  shows  the  same  signal 
with  one  digit  separation,  and  the  signal  arrival  times  are  still  resolvable. 

The  next  figure,  Figure  5.4,  demonstrates  the  processing  gain  from  coherent 
averaging.  The  data  is  the  same  as  that  in  Figure  5.1,  but  with  each  segment 
output  averaged  over  three  frames,  and  a  larger  shift  in  arrival  time.  Note  the 
increased  SNR  as  compared  with  Figure  5.1. 

Figures  5.5  and  5.6  show  the  detection  of  two  Doppler  shifted  signals.  The 
first  is  for  an  opening  situation  with  -1.4  knots  of  Doppler  and  -12  dB  SNR.  The 
second  is  a  closing  situation  with  3.4  knots  Doppler  and  a  SNR  of  -15  dB. 

Graphs  of  all  Doppler  bins  searched  are  include  in  Appendix  A.  Appendix  B 
gives  a  brief  description  of  the  program  SEQREM. 
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Figure  5.1:  Zero  Doppler  signal.  N=255,  Q=5,  fc=57  Hz,  0  dB  SNR. 

ARRIVAL  TIME  PLOT 


Figure  5.2:  Zero  Doppler  signals,  27  digits  apart.  N=255,  Q=5, 

fc=57  Hz,  -15  dB  SNR. 
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ARRIVAL  TIME  PLOT 
0.00  KNOTS  DOPPLER 


Figure  5.3:  Zero  Doppler  signals,  1  digit  apart.  N=255,  Q=5,  fc=57 

Hz,  0  dB  SNR. 

ARRIVAL  TIME  PLOT 
0.00  KNOTS  DOPPLER 


Figure  5.4:  Zero  Doppler  signal.  N=255,  Q=5,  fc=57  Hz,  -15  dB 
SNR  coherently  averaged  over  3  frames. 
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Figure  5.5:  -1.4  knot  Doppler  signal.  N=255,  Q=5,  fc=57  Hz,  -12  dB 

SNR. 


ARRIVAL  TIME  PLOT 
3.50  KNOTS  DOPPLER 


Figure  5.6:  3.4  knot  Doppler  signal.  N=255,  Q=5,  fc=57  Hz,  -15  dB 


B.  CONCLUSIONS 

The  results  of  the  previous  section  show  that  the  software  developed  for  this 
thesis  operates  as  desired.  Sequence  detection  is  performed  as  predicted,  via 
software  compensation  of  sampling  frequency,  for  any  reasonable  Doppler  using 
linear  interpolation  techniques.  Linear  interpolation  may  not  be  the  optimum 
interpolation  method,  but  it  has  been  demonstrated  that  it  does  work  and  is 
simple  and  faster  than  more  complex  techniques.  This  software  will  enable  the 
Naval  Postgraduate  School  (NPS)  to  participate  in  the  Heard  Island  Experiment 
and  will  provide  an  important  processing  capability  for  any  future  experiment, 
and  most  importantly,  a  Doppler  capability  that  was  not  previously  possessed  by 
the  Underwater  Acoustics  Systems  Laboratory  of  the  Electrical  and  Computer 
Engineering  Department  at  NPS. 

C.  ADDITIONAL  WORK 

Despite  all  of  the  advantages  of  this  software,  further  work  is  necessary. 
The  program  as  it  stands,  operates  close  to  real  time  on  individual  segments, 
given  the  inherently  low  data  rates  used  in  acoustic  tomography,  but  resides  on  a 
ULTRIX-32  V2.0  operating  system  with  a  VAX  11/785  processor.  This  is  a 
time  sharing  system  and  is  not  designed  for  real  time  processing.  Because  of 
this,  long  data  sets  end  up  with  a  reduced  priority  as  time  progresses,  causing 
slower  execution  from  the  user's  perspective. 

This  problem  will  be  remedied  with  the  delivery  of  a  Concurrent  Computer 
Corporation  VME  Native-Mode  Data  Acquisition  System  in  November  of  1990. 
The  intent  is  to  transfer  this  software  from  the  ULTRIX-32  system  to  this 
machine.  The  VME  supports  multiple  processors  and  can  process  14  MIPs  at  33 


39 


MHz.  It  uses  a  memory  cached  system  of  up  to  128  MB.  It  can  also  support  up 
to  7  separate  graphics  heads.  This  system  is  capable  of  real  time  operation  with  a 
UNIX  type  operating  system  and  will  be  ideal  for  use  with  this  programming 
package.  Given  this  system,  the  work  that  needs  to  be  done  includes: 

1.  Transfer  of  the  existing  program  and  related  functions  to  the  VME 
system.  This  includes  the  SEQREM  program,  initialization 
functions,  FHT  function,  and  magnitude  computation  routines. 

2 .  Modification  of  the  programs  to  incorporate  real  time  processing  in 
a  parallel  vice  serial  structure.  Data  input  to  the  system  is  first 
read,  processed,  stored  to  file,  and  then  the  next  segment  is  read 
and  so  on.  A  parallel  scheme  is  needed  to  allow  reading  of  the  next 
data  set  while  computation  is  progressing  on  the  current  while  the 
results  of  the  last  set  are  being  displayed. 

3 .  Design  of  an  I/O  system  to  handle  data  input  and  output  channel(s). 
Without  the  VME  available,  it  was  necessary  to  assume  input  data  is 
maintained  in  text  file.  Output  is  also  to  a  file  with  sequential 
Doppler  searches  appended  to  the  same  data  file.  The  software  also 
assumes  single  channel  input.  It  would  be  desirable  to  handle 
multiple  channels. 

4 .  Development  of  a  display  software/system  to  present  the  data  in  an 
convenient  form,  such  as  a  waterfall  type  display.  Currently, 
plotting  of  data  is  performed  off-line,  after  all  computations  are 
completed.  There  is  little  flexibility  for  plotting  various  size  data 
sets  and  real  time  display  is  not  possible. 

5 .  Altering  the  Doppler  bin  search  from  a  sequential  to  a  parallel 
operation.  Currently  in  the  data  processing  section,  Doppler 
processing  is  performed  in  increments  via  a  standard  "for"  loop. 
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The  results  of  Doppler  processing  for  the  three  possible  cases  are  given  in 
the  following  figures.  Figure  A.l  shows  a  zero  Doppler  signal.  No  scan  is 
made,  since  no  motion  is  expected.  Figure  A.2  is  a  closing  Doppler  situation. 
Doppler  search  is  over  a  +  or  -  3.5  knot  range,  in  0.5  knot  steps.  Figure  A.3  is 
the  opening  case  with  a  search  over  +  or  -  3.0  knots,  also  in  0.5  knot  increments. 
In  each  figure  the  following  parameters  apply:  primitive  polynomial  =  537s,  N  = 
255,  carrier  =  57  Hz,  \|/  =  45  degrees,  initial  state  =  lg,  Q  =  5.  White  Gaussian 
noise  is  assumed.  Plots  are  displayed  as  arrival  time  versus  time.  This  data 
simulates  a  signal  to  be  transmitted  in  the  Heard  Island  Experiment. 
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Figure  A.l:  Zero  Doppler  signal.  SNR=0  dB. 
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ARRIVAL  TiME  PLOT 
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Figure  A.2  :  3.4  knot  Doppler  signal,  0.5  knot  steps.  SNR 


ARRIVAL  TIME  PLOT 
2.00  KNOTS  DOPPLER 


(d) 

Figure  A.2  :  (cont.) 


44 


(h) 

Figure  A.2  :  (cont.) 
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ARRIVAL  TIME  PLOT 
-1.00  KNOTS  DOPPLER 


(j) 

Figure  A.2  :  (cont.) 
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(n) 

Figure  A.2  :  (cont.) 
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(C) 

S..3  :  (cont.) 


(e) 

Figure  A. 3  :  (cont.) 
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Figure  A.3  :  (cont.) 
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The  following  programs  make  up  the  SEQREM  program  or  are  associated  utility 
programs.  Since  sufficient  documentation  is  included  with  each  program,  only  a  very 
brief  description  is  provided  for  each.  This  software  may  be  obtained  on  floppy  disk 
by  contacting  Professor  James  H.  Miller  at  (408)  646-2384  or  by  writing  to  his  address 
in  the  initial  distribution  list  on  page  98. 

A.  MAIN  PROGRAM 
1.  MACROFILE 

This  is  a  a  short  file  that  holds  definition  macros  and  global  variables. 


/*  MACROFILE  includes  the  necessary  libraries  and  declares  global  constants 
and  global  variables  for  general  use. 

Variables : 

scram  -  pointer  to  indices  to  scramble  input  data  vector  prior  to  FHT. 
unscram  -  pointer  to  restore  data  vector  after  FHT  performed, 
law  -  polynomial  law  used  to  generate  M-sequence. 
degree  -  the  degree  of  law. 

initial  -  initial  register  load  for  generation  of  the  sequence.  */ 

Sinclude  <stdio.h> 
linclude  <math.h> 

Sinclude  <malloc.h> 

Sdefine  BTRMAX  11 
Sdefine  PI  3.1415926536 
Sdefine  FALSE  (unsigned) 0 
Sdefine  TRUE  (unsigned) ! FALSE 

unsigned  *scram,  *unscram,  law,  degree,  initial; 
char  *malloc(); 


2.  SEQREM 

This  program  is  the  root  of  all  the  processing  functions  for  sequence  removal. 
It  uses  standard  file  I/O  and  may  be  used  with  any  standard  input  numerical  format. 
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The  program  call  is  SEQREM,  and  the  user  is  prompted  for  the  following 
information  before  processing  may  begin  (Note:  some  parameters  are  requested  from 
within  initialization  functions,  however,  this  is  transparent  to  the  user): 

1.  Data  file  containing  formatted  filter  coefficients. 

2.  Primitive  polynomial  defining  the  m-sequence. 

3.  Initial  register  load. 

4.  Transmission  direction,  i.e.  forward  or  reverse  code. 

5.  Carrier  frequency. 

6.  Sampling  frequency  (checks  for  four  times  carrier). 

7.  Phase  modulation  angle. 

8.  Cycles  per  digit,  Q. 

9.  Desired  time  decimation,  if  any. 

10  Coherent  averaging  desired,  if  any. 

11.  Expected  Doppler  range  (knots). 

12.  Doppler  search  step  size  (knots). 

13.  Input  data  file. 

14.  Output  data  file. 

Assumptions: 

1.  Sampling  at  four  times  carrier. 

2.  One  period  of  m-sequence  processed  at  a  time. 

3.  M-sequence  is  taken  from  the  least  significant  register. 

4.  Decimation  ratio  is  a  rational. 

/*  SEQ_REM  performs  sequence  removal  of  a  phase  encoded  M-Sequence.  It  is  designed  to  work 
for  any  maximal  length  sequence  as  defined  by  a  primitive  polynomial.  The  polynomial, 
carrier  frequency,  sampling  frequency,  etc.  are  input  by  the  user.  All  filtering  and 
demodulation  is  performed,  but  filter  coefficients  must  be  stored  in  an  external  file  in  the 
format  described  by  the  function  GET_FLT_COEF() .  For  ease  of  use,  these  may  be  generated  off 
line  by  the  MATLAB  program  SET_FILTER.M. 


57 


M-Sequences  my  be  transmitted  in  the  forward  or  reverse  direction. 

?  -  is  used  as  a  wild  card  for  brevity  and  may  indicate  l,2,a,b  etc. 
in  the  following  variable  and  file  descriptions. 

Variables: 

c,  j,i,k,n,  count  -  integer  counters  used  in  various  loops. 
dir,yesno  -  variables  used  for  user  inqueries. 

num_pts  -  number  of  points  processed  per  digit,  before  and  after 
decimation . 

num_coh_avg  -  number  of  frames  to  be  coherently  averaged. 

vect_len  -  number  of  points  in  a  sequence  before  decimation  in  time. 

vect_len_inter  -  number  of  points  need  to  be  read  for  interpolation 
of  one  segment  length  at  the  new  sampling  freq. 

step  -  step  size  determined  for  decimation  in  time. 

cycles  -  number  of  carrier  cycles  per  digit,  must  be  integer  number. 

first_run  -  flag  to  determine  if  a  segment  is  the  first  one  to  be 
processed  in  the  current  doppler  bin. 

end  -  flag  to  determine  when  end  of  file  has  been  reached. 

seq_len  -  number  of  digits  in  a  H-sequence. 

highpass_?, lowpass?_?  -  pointer  to  storage  for  filter  coefficients. 

indatal, indata2  -  array  for  data  storage  and  processing  storage  before 
decimation  in  time.  Indatal  is  initially  input 
and  then  is  imaginary  part  and  indata2  is  real  part 
(after  demodulation) . 

h_data_im, h_data_re  -  arrays  for  real  and  imaginary  vectors  for  storing 
scrambled  data  for  FHT  processing. 

re__data,  im_data  -  arrays  for  real  and  imaginary  data  storage  after  data 
unscrambling. 

avg_re,  avg_im  -  arrays  for  real  and  imaginary  data  storage  when 
coherent  averaging  is  performed. 

mag,  phase  -  array  for  magnitude  and  phase  of  processed  data. 

fc, fs  -  pointers  for  carrier  frequency  and  sampling  frequency. 

dclvl_im,dclvl_re,pdstal_re,pdstal_im  -  used  to  remove  correlation 

dc  bias. 

ifact,rfact,ph_ang  -  factors  and  phase  angle  used  to  compute 
correction  to  dc  bias  based  on  phase  modulation 
angle. 

den  -  tmeporary  storage  used  in  above  calculations  of  dc  bias. 
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snd_vel  -  input  average  sound  velocity  to  be  used. 

Ts,Ta,Ti  -  sampling  period,  sequence  period,  and  interpolated  sampling 
period  in  doppler  interpolation. 

time_s,time_i  -  current  relative  sample  time  and  interpolation  time, 
which  are  used  to  compute  h. 

h  -  distance  from  time  a  to  sample  to  be  interpolated  where 
time  a  <  h  time  b. 

ptrl,ptr2  -  used  to  save  position  in  array  space  when  performing 
doppler  search. 

doppler  -  current  doppler  shift  being  processed  (m/s) . 

doppler_start,doppler_end  -  beginning  and  ending  of  doppler  search 
interval  (m/s)  . 

doppler_step,doppler_bin  -  doppler  step  size  used  in  scan,  and  doppler 
converted  to  radians  for  demodulation. 


Functions: 

init_had()  -  gets  relevant  parameters  for  computing  scambling  and 
unscrambling  indices. 

fwd_had() ,  rev_had()  -  compute  indices  for  foward  and  reverse  trans¬ 
mitted  sequences.  Both  return  polynomial 
degree  from  input. 

lowpass?  0 , highpass ()  -  identical  functions  for  performing  filtering 
operations. 

get_flt  coef()  -  retrieves  filter  coefficient  data  from  user  provided 
"file. 

demodulateO  -  performs  demodulation  for  input  signal  to  dc.  Has  an 
argument  dopplerjbin  for  doppler  frequency  shifts. 

hadamardO  -  performs  the  FfiT  on  scrambled  data.  Identical  to  FFT 
with  multiplication  factors  omitted. 

mag_phase()  -  computes  the  magnitude  and  phase  of  the  processed  data.*/ 


f include  "macrofile. h" 

main  () 

{ 

int  c,  j,k,n,  i, dir,  cycles, num_pts,  step, vect_len,vect_len_intor, count, num_coh_avg; 
unsigned  init_had(),  rev_had() ,sec_len,  first_run,  end; 

void  get_flt_coef  () , hipass () , lowpassl  () , demodulateO , hadamardO  ,mag_phase{)  ; 
void  lowpass2 0 , lowpass3() ; 

float  *hipass_a, *hipass_b, *lowpass_a, *lowpass_b, *lowpassl_a, "lowpassljb; 
float  *indata2,  *h_data_re,  *h_data_im,  »re_dara,  *ira_data,  *fs,  »fc; 
float  *roag,  'phase, *avg_re, *avg_in,dclvl_re,dclvi_im,pdscal_re,pdstal_im; 
float  *data_i_ro,  *data_i_ira,  rfact,  ifact,  den,  ph_ang,  doppler,snd_vel; 
float  doppler_start,time_s,  time_i,  h,  Ts,  Ta,  Ti,  *tcmp_ptri,  *tenp_ptr2; 
float  doppler_end,doppler_bin,doppler_stcp,  »indatal,  *nevf()  ; 
int  yesno; 


59 


char  infile(21],  outfilel21); 
FILE  *fpl, *fp2; 


/*  Initialize  filter  coefficients  for  a  low  and  high  pass  filter.  */ 

/*  These  values  are  computed  and  stored  in  a  file  off  line.  */ 

count  =0; 

/*  Allocate  memory  Tor  fc,  fs,  and  filter  coefficients  to  be  used.  10th 
order  BUTTERWC"' bandpass,  and  5th  order  CHEBYCHEV  lowpass  filters 
must  be  used.  Filter  coefficients  are  retrieved  and  filters  are 
initialized.  Memory  allocated  to  filter  coefficients  is  freed  after 
initialization.  */ 

fc  =  new(l)  ; 
f  s  =  new(l)  ; 
hipass_a  **  new(BTRMAX); 
hipass_b  new(BTRMAX) ; 
lowpass_a  »  new(BTRMAX); 
lowpass_b  °  new (BTPMAX) ; 
lowpassl_a  °  new(BTRMAX); 
lowpassl_b  »  new (BTRMAX) ; 

get_flt_coef (hipass_a,hipass_b, lowpass_a, lowpass_b, lowpassl_a,  lowpassl_b) ; 

hipass (hipass_b, hipass_a, BTRMAX) ; 

lowpassl.(lowpass_b,  lowpass_a,  BTRMAX-5) ; 

lowpass2  (lowpassjo,  lowpass_a,  BTRMAX-5) ; 

lowpass3 (lowpassl_b,  lowpassl_a, BTRMAX) ; 

free(hipass_a) ; 

free (hipass_b) ; 

free (lowpass_a) ; 

free (lowpass_b) ; 

free  (lowpass,l_a) ; 

free (lowpass l_b) ; 

/*  Set  sound  velocity  to  desired  value,  for  given  conditions.  */ 

print; ("\nEnter  the  Average  Sound  Velocity  (muters/sec. ) . ") ; 
printf("\n\n  Sound  Velocity:  "); 

scant  ("%f ",  «snd_vel) ; 

/*  Initialize  parameters  for  M-Sequence,  and  set  transmit  direction. 

S^-t  scrambling  and  unscrambling  indices  accordingly.  */ 

initial  •*  init_had(); 

printf("\nls  the  M-Sequence  transmitted  In  the  toward  or  reverse\n,:)  ; 
printf ("direction?  \n"); 

printf("\n  (forward=0/reverse=l) :  "); 

scanf  ("%d",  sdir) ; 
if  (dir  ==  0) 

degree  “  fwd_had(law, initial, scram, unscram); 

else 

degree=  rev_had(law,  initial, scram, unscram) ; 
seq_len  -  (l«degree)-l; 


/*  Compute  offset  to  remove  DC  bias  from  correlation.  */ 

print; ("\n3ntcr  the  Phase  Modulation  Angle  used  to  encode  the  signal. \n  ") 
printf ("\n  Phase  Angle:  ") ? 

scanf {"%£", Sph_ang) ; 
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ph_ang  =  ph_ang*PI/180; 

den  >=  seq_len*seq_len  +  tan  (ph_ang)  *tan  (ph_ang) ; 

ifact  =  (seq_len+l) *tan (ph_ang) /den; 

rfact  =  - (seq_len-tan (ph_ang) *tan (ph_ang) ) /den; 

/*  Input  of  transmission  parameters  and  demodulator  is  initialized.  */ 

printf  ("\nEnter  the  Carrier  Frequency  usedAn  ") ; 
printf ("\n  Fc:  "); 

scanf ("%f", fifc(O) ) ; 

printf ("\nEnter  the  Sampling  Frequency.  Sampling  Frequency  must  be  an"); 

printf  ("\nfour  times  the  Carrier  Frequency  An")  ; 

printf ("\n  Fs:  ") ; 

scanf  (:,%f",  Sfs  (0) ) ; 

if  (  *fc/(*fs)  !=  0.25) 

{ 

printf ("\nSampling  Frequency  will  not  work  with  this  program.  Exiting! \n") ; 
exit (1) ; 

) 

printf ("\nEnter  the  number  of  carrier  cycles  per  digit.  An  "); 
printf ("\ninteger  number  of  cycles  must  be  usedAn") ; 
printf ("\n  Q:  "); 

scanf ("id", scycles) ; 

/*  Compute  sequence  length  in  seconds  (Ta)  and  sampling  period.  The  number 
of  points  generated  per  digit  is  computed  and  provided  as  an  aid  in 
determining  the  decimation  in  time  to  be  used,  if  desired.  The  total 
number  of  points  per  sequence  period  is  computed  (vect_len) . */ 

num_pts  =  (*fs/ (*fc) ) *cycles; 
vect_len  »  num_pts*seq_len; 

Ta  =  (float) (seq_len  *  cycles  /  *fc);/*  computes  the  sequence  length  in  sec.*/ 

Ts  “  l/(*fs);  /*  compute  sampling  period.  */ 

printf ("\nThere  are  %d  pts  per  digit.  For  Processing  Savings  ",num_pts); 
printf ("\ndecimation  in  time  is  performed.  The  default  used  is  one  point"); 
pri  ltf ("\nfor  each  cycle  in  the  digit.  In  this  case  %d  points  are", cycles) ; 
printf ("\nprocessed  per  digit,  (i.e.  every  %d  th  point .)", (num_pts/cycles) ) ; 
printf("\n\n  Use  default? (yes=l/no=0) :  "); 

scanf ("%d", Syesno) ; 
if  (yesno  ==  1  ) 

step  =  num_pts/cycles; 
else 
{ 

printf ("\nEnter  the  desired  Decimation.  (e.g.  for  every  point  enter  l,\n"); 
printf ("for  every  other  point  enter  2,  for  every  third  point  enter  3,  etc"); 
printf ("\nDecimation  must  be  evenly  divisible  into  the  pts.  per  digit."); 
printf ("\n (NOTE:  Decimation  must  be  less  than  %d  tc  ensure  \n", num_pts) ; 
printf  ("  at  least  one  point  per  digit  is  used!)\n"); 

printf ("\n  Decimation:  ") ; 

scanf ("%d", Sstep) ; 
if  (  num_pts%step  !=  0) 

( 

printf ("\nlnvalid  decimation.  Try  again."); 

printf ("\n  Decimation:  "); 

scanf ("%d", Sstep) ; 

if  (  num_pts%step  !=  0) 

printf ("\nError .  Bye!\n"); 

> 

if  (step  >  num_pts) 

{ 

printf ("XnDecimation  chosen  is  too  large,  use  a  smaller  number. \n"); 
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printf("\n  Decimation:  "); 

scanf ("%d", &  step) ; 
if  (step  >  num_pts) 

ptintf ("\nSorry !  Still  won't  work,  aborting\n") ; 

) 

} 

/*  Determine  if  coherent  averaging  of  sequences  is  desired  for  improving 
processing  gain.  */ 

printf ("\nEnter  the  number  of  Sequence  Frames  to  be  coherently  averaged."); 
printf  ("\nNumber  Frames  »  1  implies  no  coherent  averaging  desired. \n\n") ; 
printf ("  Number  Frames:  ")  ; 

scanf ("%cl", snum_coh_avg) ; 
if  (num_coh  avg  <=  0) 

( 

printf ("\nlnvalid  Number  of  Frames.  Must  Use  Positive  Integer."); 
printf ("\n\n  Number  Frames:  ") ; 

scanf ("%d", tnum_coh_avg) ; 
if  (num_coh_avg  <“  0) 

< 

printf ("\nlnvalid  Number  of  Frames.  Aborting ! \n") ; 
exit (1) ; 

) 

) 

/*  Query  for  Doppler  range  to  be  searched.  An  input  of  zero  will  cause 
the  resampling  for  doppler  steps  to  be  skipped.  Doppler  is  coverted 
from  knots  to  meters/sec  for  computation.  Sequence  frequency  resolution 
is  computed  to  assist  in  determining  doppler  bins  to  be  searched.  */ 

printf ("\nEnter  the  doppler  range  to  be  searched  (knots)."); 
printf ("\n\n  Doppler  (+/-):  "); 

scanf ("%f", Sdoppler) ; 

printf ("\nOne  knot  corresponds  to  a  %1.5f  Hz  shift", *fc-*fc*(l. 0-0. 5/snd_vel) ) ; 
printf ("  in  frequnecy . \nEnter  the  Doppler  increment  to  be  searched."); 
printf  ("\n\n  Step  (knots):  "); 

scanf ("%f", Sdoppler_step) ; 
if  (doppler_step  ““  0.0) 
doppler_step  =  1.0; 

doppler  /=  2.0;  /*  convert  knots  to  meters/second  */ 

doppler_step  /=  2.0; 
doppler_start  -  -doppler; 
doppler_end  =  doppler; 

/*  Query  for  input  and  output  files.  */ 

printf ("\nEnter  Input  File  Name  (20  Characters  Maximum) . \n") ; 
printf ("\n  File  Name:  "); 

scanf ("%s", inf ile) ; 

printf ("\nEnter  Output  File  Name  (20  Characters  Maximum) .\n") ; 

printf  ("\n  File  Name:  ”); 

scanf ("%s", outfile)  ; 

if  ((fpl  =  f open (infile, "r") ) ==NULL) 

{ 

printf ("\nUnable  to  Open  Input  File.  Aborting ! \n") ; 
exit  (1)  ; 

) 

if  ((fp2  “  fopen (outfile, "w") ) =“NULL) 

{ 

printf ("\nUnable  to  Open  Output  File.  Aborting ! \n") ; 
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exit  (1) ; 
} 


/*  After  Decimation  there  are  less  data  pts  =>  comput  new  num_pts.  */ 
num_pts  /**  step; 

/*  Set  data  length  for  input  to  maximum  needed  for  highest  doppler, 

which  corresponds  to  the  maximum  number  of  points  used  in  interpolation.*/ 

Ti  =  Ts/(1  +  doppler_end/snd_vel)  ; 
vect_len_inter  =  Ti*vect_len/Ts  +  1; 

/*  Allocate  memory  as  required.  */ 

indatal  =  new(vect_len_inter) ; 
indata2  =  new (vect_len_inter) ; 
im_data  =  new(num_pts*seq_len) ; 
re_data  «*  new(num_pts*seq_len) ; 
mag  *■  new(num_pts*seq_len) ; 
phase  =  new (num_pts*seq_len)  ; 
avg_re  =  new(num_pts*seq_len) ; 
avg_im  =  new(num_pts*seq_len) ; 
h_data_im  =  new(seq_len+l) ; 
h_data_re  =  new(seq_len+l) ; 
data_i_re  =  new(vect_len) ; 
data_i_im  =  new(vect_len) ; 

/*  DATA  PROCESSING  SECTION. 

Doppler  processing  is  performed  automatically  according  to  the  step 
size  specified  by  the  user.  Data  is  processed  in  lengths  corresponding 
one  segment  at  a  time.  Sufficient  points  for  complete  interpolation 
of  one  doppler  shifted  sequence  are  picked  off  the  input  file. 

During  doppler  shift  interpolation  the  input  end  point  is  saved  for 
processing  in  the  next  segment  to  insure  segment  continuity.  */ 

for  (doppler  =  doppler_start; doppler  <=  doppler_end;  doppler+=doppler  step) 

{ 

/*  Initialize  demodulation  routine.  */ 
demodulate (fc, fs, doppler, vect_len) ; 

/*  Determine  the  sampling  period  for  the  new  sampling  frequency  used  in 
interpolation,  as  a  function  of  doppler.  */ 

Ti  *»  Ts/(1  +  doppler/snd_vel) ; 

/*  If  zero  doppler  case  is  being  performed  endpoint  bookkeeping  is  not 
used  and  data  length  is  the  same  as  input  data  (vect_len) .  */ 

if  (doppler  !**  0  ) 

< 

fscanf (fpl, "%f\n",  Sindatal (0) )  ; 
vect_len_inter  =  Ti*vect_len/Ts; 
indatal++; 

*indata2++  =  0; 

) 

else 

vect_len_inter  =  vect_len; 

/*  Initialize  flags  and  interpolation  parameters.  t'irst_run  indicates 
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first  segment  being  processed,  and  end  indicates  EOF  reached.  Time_s 
indicates  current  relative  sampling  time  and  time_i  indicates  current 
relative  interpolation  time.  */ 

time_s  =  0.0; 
time_i  =  Ti; 
end  =  FALSE; 
first_run  «=  TRUE; 

dopplerjbin  =  PI  *  doppler/ (2*snd_vel) ; 

/*  Ta  is  sequence  length  in  seconds.  Header  is  printed  to  separate 
doppler  scans.  */ 

Ta  =  cycles  *  seq_len  /(*fc  *  (1  +  doppler/snd_vel) ) ; 
fprintf (fp2, "START\n") ; 

fprintf (fp2, " 1 %6.2f  KNOTS  DOPPLER' \n", 2*doppler) ; 

fprintf (fp2,"' SEQUENCE  PERIOD  (%7.4f  SEC. ) ■\n",Ta) ; 

fprintf (fp2, "%5d  %10.7f\n", seq_len*num_pts, Ta/ (seq_len*num_pts-l) ) ; 

/*  Process  entire  input  file  for  current  doppler  search.  */ 

while  (lend) 

{ 

for  ( j=0;_j<vect_len_inter;  j++) 
if  ( (c=fgetc(fpl)7==E0F) 

{ 

rewind (fpl) ; 

demodulate (fc,fs, doppler, 0) ;  /*  clears  demodulator  for  next  run.*/ 
end  =  TRUE; 
break; 

) 

else 

{ 

ungetc(c, fpl) ; 

f scanf (fpl, "%f\n", sindatal 1 j) ) ; 

) 

/*  Test  for  EOF  reached  before  complete  segment  read.  Prevents  processing 
partial  segments.  */ 

if  (lend) 

{ 

/*  Test  for  first  segment.  If  it  is  the  first  data  point  has  not  been 
filtered,  adjusts  bookkeeping  of  endpoint.  */ 

if  ( (doppler 1=0.0)  &&  first_run) 

{ 

— indatal; 

— indata2; 

vect_len_inter  +=  1; 

) 

/*  Lowpassing  and  then  highpassing  reduces  to  a  BP  equivelant  and  avoids 
complex  numbers  at  this  point.  Caution:  the  filter  phase  must  be  linear 
in  the  pass  region.  */ 

hipass (indatal, indatal, vect_len_inter) ; 

lowpass3 (indatal, indatal, vect_len_inter) ; 

demodulate ( indatal, indata2, doppler_bin, vect_len_inter ) ; 

lowpassl (indatal, indatal, vect_len_inter)  ; 

lowpass2 (indata2,indata2,vect_len_inter)  ; 
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/*  This  loop  perforins  a  linear  interpolation  of  the  demodulates.  Maintains 
bookkeeping  of  endpoint,  and  first  segment.  Interpolates  vect_len  points 
for  hadamard  processing.  Skips  this  if  zero  doppler  case.  */ 

if  (doppler  !=  0.0) 

{ 

i=0; 

if  (first  run) 

( 

vect_len_inter  -=  1; 
first  run  =  FALSE; 

) 

else 

{ 

— indatal; 

— indata2; 

) 

for  (j=0;  j<vect_len;  j++) 

{ 

/*  First  case  is  down  shift  of  doppler.  */ 

if  (time_i  >  time_s+Ts) 

( 

time_s  +*=  Ts; 
i++; 

h  “  time_i  -  time_s; 

data_i_im(j)  =indatal(i]  +  (h/Ts) * (indatal [i+l)-indatal (i) ) ; 
data_i_re(j]  =indata2(i)  +  (h/Ts) * (indata2 [ i+1 ) -indata2 [ i) ) ; 

/*  Reset  sampling  time  and  interpolation  time.  The  relative  time  is  important 
only.  */ 


time_s  =  0.0; 
time_i  »  h  +  Ti; 
if  (time_i  >  2*Ts) 

( 

time_s  +=  Ts; 

i++; 

) 


/*  Second  case  is  upshift  in  doppler.  */ 

else 

< 

h  =  time_i  -  time_s; 

data_i_im{j]  =  indatal(i)  +  (h/Ts) * (indatal (i+1 ) -indatal (i) ) ; 
data_i_re(j)  =indata2(i]  +  (h/Ts) * (indata2 [i+1 ) -indata2 [i] ) ; 
time_i  +=  Ti; 
if  (time  i  >  time_s  +  Ts) 

( 

time_s  +**  Ts; 

i++; 

time_i  -=  time_s; 
time_s  =  0.0; 

) 

} 

) 
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/*  Save  indatal  and  indata2  working  space.  */ 

temp_ptrl  =  indatal; 
temp_ptr2  =  indata2; 

*temp_ptrl++  =  indatal (i+1] ; 

*temp_ptr2++  =  indata2 (i+1) ; 

/*  Set  pointers  to  interpolated  values  */ 

indatal  =  data_i_im; 
indata2  =  data_i_re; 

} 

/*  Decimate  in  time  according  to  step.  And  compute  sequence  dc  level.  */ 
n=0; 

dclvl_re  =  0.0; 
dclvl_im  =  0.0; 
for  ( j=0; j<vect_len; j+=step) 

{ 

im_data(n)  =  indatal (jj; 
re_data(n++]  =  indata2(j); 

}  “ 


/*  Restore  indatal  and  indata2  working  space,  when  processing  for  doppler.  */ 


if  (doppler 
{ 

indatal  = 
indata2  = 
} 


!=  0.0) 

temp_ptrl; 

temp_ptr2; 


/*  Scramble  data  vector  and  process  data  according  to  interleave  via  FHT 
(HADAMARD)  and  unscramble.  Real  and  imaginary  parks  are  processed 
separately.  */ 

for  ( j“0; j<num_pts; j++) 

( 

h_data_im ( 0 )  =  0 ; 
h_data_re(0)  =  0; 
n=0; 

for  (k=0;  k< (seq_len*num_pts) ;k+»num_pts) 

{ 

h_data_im (scram (n ] )  «*  im_data(k+ j] ; 
h_data_re ( scram ( n++ ) )  =  re  data(k+j); 

) 

hadamard (degree, h_data_ira) ; 
hadamard (degree, h_data_re) ; 
n=0; 

dclvl_re+=h_data_re (0] ;  /*  Save  pedestal  information.  */ 

dclvl_im+=h_data_im(0] ; 

for  (k=0;k<(seq  len*num_pts) ;k+=num_pts) 

{ 

im_data[k+j)  =  h_data_im(unscram(n)) ; 
re_data(k+j)  =  h_data_re (unscram(n++) ) ; 

) 

) 


/*  Compute  pedestal  corrections  for  real  and  imaginary  parts.  */ 

pdstal_re  =  (dclvl_re/num_pts) *rfact  -  (dclvl_im/num_pts) *ifact; 
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pdstal_im  =  (dclvl_re/num_pts) *ifact  +  (dclvl_im/numjpts) *rfact; 

/*  Case  1  is  if  Coherent  averaging  is  performed.  Else,  no  coherent  averaging 
is  performed.  Each  segment  is  output.  */ 

if  (num_coh_avg  >  1) 

{ 

/*  Make  DC  level  correction.  */ 

for (k=0;k<seq_len*num_pts;k++) 

( 

avg_re[k)+=  (re_data[k)-pdstal_re) ; 
avg_im(k]+=  (im_data[k)-pdstal_im) ; 

} 

count++; 

if  (count==num_coh_avg) 

{ 

count  =0; 

for  (k“0;k<seq_lon*num_pts;k++) 

( 

avg_re ( k ) /=num_coh_avg; 
avg_im ( k ] / =num_coh_avg ; 

) 

/*  Compute  magnitude  and  phase  and  print  results.  */ 

mag_phase(avg_re,avg_im, mag, phase, (num_pts*seq_len) ) ; 
for  ( j°0; j<num_pts*seq_len; j++) 

{ 

fprintf (fp2, "%8. If  %8 . If \n", mag [ j) , phase ( j) ) ; 
avg_re[j]  =  0.0; 
avg  im[ j  J  =  0.0; 

) 

) 

) 


/*  Case  2.  No  coherent  averaging.  Compute  magnitude  and  phase  and 
print  results.  */ 

else 

{ 

for  (k=0; k<seq_len*num_pts; k++) 

( 

re_data[k)-=  pdstal_re; 
im_data[kJ-=  pdstal  im; 

} 

mag_phase(re_data,im_data, mag, phase, (num_pts*seq_len) ) ; 
for  ( j“0; j<num_pts"seq_len; j++) 

fprintf {fp2, "%8.1f  %8 .lf\n",mag l j) , phase ( j] ) ; 

) 

) 

} 

) 

/*  Finished  close  out  variables  and  files.  */ 

f close (fpl) ; 
f close (fp2) ; 
free(fc)  ; 
free (fs) ; 
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free (indatal) ; 
free (indata2) ; 
free (im_data) ; 
free (re_data) ; 
free (mag) ; 
free (phase) ; 
free (avg_re) ; 
free  (avg_im)  ; 
free (h_data_im) ; 
free (h_data_re) ; 
free (data_i_re) ; 
free (data_i_im)  ; 

) 

/*  NEW  is  a  short  function  to  allocate  memory  for  floating  point  vectors.  */ 

float  *new(size) 
int  size; 

( 

float  ‘newdata; 

if  ((  newdata  =  (float  *)malloc(size*sizeof (float) ) )==NULL) 

( 

printf ("Cannot  Allocate  Storage !\n"); 
exit  (1)  ; 

) 

return (newdata) ; 

) 


B.  INITIALIZATION  PROGRAMS 
1.  INIT  HAD 

This  program  prompts  for  the  primitive  polynomial  initial  register  load  and 
allocates  memory  space  to  the  variables  scram  and  unscram.  Returns  initial  register 
load,  initial. 

/*  INIT_HAD  is  a  program  to  initialize  memory  allocation  for  the  scrambling 
and  unscrambling  arrays  to  be  used  with  Fast  Hadamard  Transform.  It 
also  returns  the  initial  register  load  used  with  the  Shift  Register 
Generator. 

Variables : 

scram  -  external  array  to  hold  scrambling  indices. 

unscram  -  external  array  to  hold  unscrambling  indices. 

initial  -  initial  register  load  for  SRG. 

law  -  polynomial  law  that  defines  the  SRG  structure. 

length  -  length  of  sequence  generated  by  the  input  law.  Assumes 
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maximal  length  polynomial. 

temp  -  temporary  holding  spot  for  determining  sequence  length.*/ 


linclude  "macrofile.h" 
unsigned  init_had ( ) 

{ 

extern  unsigned  law; 
unsigned  initial,  temp,  length; 

printf ("\nEnter  polynomial  law  for  the  desired  M-Sequence.  Use"); 
printf ("\noctal  integer  representation  only."); 
printf ("\n(e.g.  D3  +  D  +  1  -  1011  binary  =>  13  Octal.) \n") ; 
printf <"\n  LAW:  "); 

scanf ("%o", Slaw) ; 
temp  “  law; 
length  =  1; 
while  (temp»=l) 
length«=l; 
length — ; 

if  ((scram  *»  (unsigned  *) malloc (length*sizeof  (unsigned) ) )  ==  NULL) 

{ 

printf ("\nCannot  Allocate  Scram  Array! !! !\n") ; 
exit  (1) ; 

) 

if  ((unscram  ■  (unsigned  *)  malloc  (length*sizeof  (unsigned) ) )  ■**=  NULL) 

( 

printf ("\nCannot  Allocate  Unscram  Array! !! !\n") ; 
exit  (1)  ; 

) 

printf ("\nEnter  the  initial  register  load  in  Decimal,  do  not  use  zero."); 
printf ("\n\n  Initial  Lead:  ") ; 

scanf ("%u", Sinitial) ; 
if  (initial  *»■  0) 

( 

printf ("\nError!  Initial  load  cannot  be  0!!\n\n"); 

printf ("\nEnter  the  initial  register  load  in  Decimal,  do  not  use  zero."); 
printf ("\n\n  Initial  Load:  "); 

scanf ("%u", sinitial) ; 
if  (initial  =«•  0) 

( 

printf ("\nError !  Aborting.\n\n") ; 
exit  (1) ; 

) 

) 

return (initial) ; 

) 


2.  FWD  HAD 

Computes  the  permutation  indices  for  use  with  the  FHT.  Data  is  transmitted  in 
the  forward  direction.  Returns  the  degree  of  the  primitive  polynomial. 

/*  FWD_HAD()  computes  the  permutation  indices  for  scrambling  and 

unscrambling  a  data  vector  generated  using  Maximal  Length  Sequences. 

These  indices  are  used  with  the  Fast  Hadamard  Transform  and  is  performed 
in  place.  Sequences  are  assumed  to  be  transmitted  in  the  foward 
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direction . 


Variables: 

law  -  Shift  Register  Generator  law  to  be  used. 

initial  -  initial  register  law. 

scram  -  array  of  indices  for  scrambling  data. 

unscram  -  restoration  indices  for  use  after  FHT. 

degree  -  gives  the  degree  of  the  polynomial  law. 

temp  -  temporary  holder  for  computing  degree  and  seq_len. 

S_law  -  law  formed  from  law  to  implement  the  S_gen  structure. 

rev_law  -  law  in  the  reverse  order  to  implement  the  L_gen  structure. 

end_bit  -  maintains  end  bit  for  around  end  feed  in  S_gen. 

L_gen  -  variable  that  acts  as  L  generator  delay  registers. 

S_gen  -  variable  that  acts  as  S  generator  delay  register. 

Reference:  The  Feedback  generators  used  for  forming  the  foward  scrambling 
and  unscrambling  indices  are  adapted  from: 

M.  Cohn  and  A  Lempel,  'On  Fast  M-Sequence  Transforms’, 

IEEE  Transactions  on  Information  Theory,  Jan.  1977.  */ 

unsigned  fwd_had(law, initial, scram, unscram) 
unsigned  law,  initial,  *scram,  *unscram; 

{ 

unsigned  degree, temp, temp2, seq_len, S_law, rev_law; 
unsigned  end_bit, L_gen,S_gen; 
int  i, j, count ; 

/*  Initialize  variables.  */ 
temp  =  law; 
degree  =0; 

S_gen  =0; 
rev_law  =  0; 
seq_len  =  1? 

/*  Computes  the  length  of  the  M-Sequence,  and  Polynomial  Degree.  */ 
while (temp>>=l) 

{ 

seq_len  <<»  1; 
degree++; 

) 

/*  Reverses  law  for  use  with  L_gen.  */ 
temp  =  law; 

for  (i=0;  i<=  degree;i++) 

{ 

rev_law  -  (tempSl) I (rev_law<<l) ; 
temp>>»l; 

} 
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/*  Set  end_b.it  for  register  end  feed  for  S_gen  logic.  Set  initial  load 
for  S_gen.  Set  law  for  use  in  S_gen  logic.*/ 
end_bit  =  seq_len — ; 

S_law  =  :(law»l)  ; 

end_bit»=l; 

temp  =  initial; 

for  (i=0;  i<seq_len-l;  i++) 

{ 

if  (temps 1) 

temp  *»  (tempAlaw)  »1; 
else 

temp  »=  1; 

if(i  >=  seq  len-degree) 

{ 

S_gen  =  S_gen|  (temps  1)  ; 

S  gen  «=1; 

} 

} 

S_gen  1=  (initialsl); 

/*  Load  L_gen  to  generate  unscrambling  indices.  */ 

L_gen  *»  (1«  (degree-1) )  ; 

/*  Compute  permutations  scram  and  unscram.  */ 
for  (i=0;  i<seq_len;  i++) 

( 

unscram(i)  =  L_gen; 

temp2  «■  0; 

temp  =  S_gen; 

for  ( j=0; j<degree; j++) 

{ 

temp2  «=1; 
temp2  1°  (temps  1); 
temp  »*1; 

) 

scram (i]  =  temp2; 
if  (L_gensl) 

L_gen  °  (L_gen/'rev_law)  >>1; 
else 

L_gen»°l; 

temp  =  (S_gen  s  S_law) ; 

/*  Count  the  number  of  l's  for  modulus  two  sum  for  end  feedback  to  S_gen.  */ 
count=0; 

for  (j=0; j<degree; j++) 

( 

if  (tempsl) 

{ 

count++; 

temp»al; 

} 

else 

temp>>=l; 

} 

if  (  (count%2)  ==  0) 

S_gen  <<»  1; 
else 

S_gen  =  (S_gen«l)|l; 

S  gen  S»  seq_len; 

)“ 

return (degree)  ; 
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) 


3.  REVHAD 

Computes  the  permutation  indices  for  use  with  the  FHT.  Data  is  transmitted  in 
the  reverse  direction.  Returns  the  degTee  of  the  primitive  polynomial  [Ref.  15]. 

/*  R£V_HAD()  computes  the  scrambling  and  unscrambling  indices  for  a 

Maximal  Length  Sequence  that  is  to  be  processed  using  a  Fast  Hadamard 
Transform  and  is  performed  in  place.  The  sequence  is  assumed  to  be 
transmitted  in  the  reverse  direction. 

law  -  The  polynomial  law  defining  the  Maximal  Length  Sequence  that  is 
used. 

initial  -  The  initial  register  load  for  the  shift  register  generator. 

scram  -  Pointer  to  the  array  to  contain  the  scrambling  indices, 
to  be  used  with  the  transformed  array.  The  transform 
is  not  done  in  place. 

unscram  -  Pointer  to  the  array  to  contain  the  unscrambling  indices. 

The  transform  is  not  done  in  place. 

degree  -  The  degree  of  the  law  is  returned. 

Reference:  This  procedure  was  adapted  from  the  article  "On  Fast 
M-Sequence  Transforms",  by  Martin  Cohn  and  Abraham  Lempel, 

IEEE  Transactions  on  Information  Theory,  Jan.  1977. 

The  program  was  modified  from  code  provided  by  Kirk 
Metzger  with  permission.  */ 

Sinclude  "macrofile.h" 

unsigned  rev_had(law, initial, scram, unscram) 
unsigned  int  law,  initial,  *scram,  *unscram; 

( 

/*  temp  Scratch  variable  for  use  in  finding  degree  etc. 

endjbit  Bit  is  set  corresponding  to  highest  bit  in  law.  (i.e.  the 
end  register  n) . 

rev_initial  Initial  register  load  for  reverse  generation. 
seq_len  Length  of  sequence  based  on  law. 

index  Loop  counter. 

ss_contents  Shift  register  contents  for  scram  array. 
ms_contents  M  sequence  register  contents  for  unscram  array.*/ 

unsigned  int  degree,  temp,  end_bit,  rev_initial,  seq_len,  index; 
register  unsigned  ss_contents,  ms_contents; 

/*  Initialize  variables.  */ 
temp  “  law; 
sec_len  ■*  1; 
rev_initial  =  0; 
degree  =  0; 

/*  Find  sequence  length  and  degree  of  polynomial.  */ 
while  (temp>>»l) 

{ 


72 


sea_len«=l; 

degree++; 

} 

/*  set  end_bit  */ 

endjbit  »  seq_len — »1; 

/*  find  reverse  generator  initial  load  */ 
for (index  =  0;  index  <  degree;  index++) 

( 

rev_initial  =  (rev_initial«l)  I  (initialsi)  ; 
initial»°l; 

} 

/*  generate  scram  and  unscram  values  using  law,  end_bit  and  values  of 
ss_contents  and  ms_contents.  */ 
ms_contents  =1; 
ss_contents  =  rev_initial; 
for  (index  =0;  index<seq_len;  index++) 

( 

*scram++  =  ss_contents; 

*unscram++  =  ms_contents; 
temp  =  ss_contentsSlaw; 
ss_contents»=l; 
do  {  if  (tempsi) 

ss_contents"'=endJbit; ) 
while  (temp>>=l) ; 
if  (ms_contents&l) 

ms_contents  =  (ms_contents,'law)  »1; 
else 

ms_contents»=l; 

) 

return (degree) ; 

) 

4.  GETFLTCO. 

Retrieves  the  filter  coefficients  to  be  used  used  in  initializing  the  filtering 
programs.  Two  tenth  order  and  one  fifth  order  filters  are  used  in  this  program. 


/*  GET_FLT_CO£F  retrieves  the  filter  coefficients  from  the  user  specified 
file  using  a  columnar  format.  The  first  column  is  the  b  coefficients 
and  the  second  is  the  a  coefficients.  The  first  11  rows  are  the  highpass 
coefficients  used  for  the  upper  end  of  the  passband  filter  and  the  next 
eleven  are  the  lowpass  coefficients  for  the  lower  end  of  the  passband. 
the  final  six  rows  are  for  the  low  pass  filter  used  after  demodulation . 

The  coefficients  file  may  be  created  by  any  program  desirea  by  the  user 
but  must  be  in  the  appropriate  format.  However,  a  MAT LAB  program, 
set_filter ,m,  does  this  automatically.  Set_filter.m  uses  a  BUTTERKOWH 
filter  (10th  order)  for  the  passband  and  a  CHEBTCHEV  filter  for  the 
lowpass  filter  (5th  order) .  •/ 


^include  <stdio.h> 
^include  <math.h> 
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void  get_flt_coe£ (hi_a, hi_b, low_a, low_b, lowl_a, lowl_b) 
float  *hi  a,  *hi_b,  *\low_a,  *low_b,  *lowl_a,  *lowl_b; 

{ 

int  i; 

FILE  *fp; 

char  filterfile(21) ; 

printf ("\nEnter  the  filename  with  the  desired  filter  coefficients."); 
printf("(20  Characters  Max)\n"); 
printf ("\n  Filename:  "); 

scar.f  ("%s",  interfile) ; 
if((fp  ■  fopen (filterfile, "r") ) ==  NULL) 
f 

printf (“\nUnable  to  open  file.  Aborting ! \n") ; 
exit  (1) ; 

) 

for  (i=0;i<=10;i++) 

fscanf (fp, "%f  %f \n", &hi_b(i) , shi_a[i) ) ; 
for  (i=0;i<=10;i++) 

fscanf (fp, "%f  %f\n",&lowl  b(i) , filowl_a (i) ) ; 
for  (i=0;  i<=5;  i++) 

fscanf (fp, "%f  %f\n", Slow_b(i) , Slow_a [ i) )  ; 
f close (fp) ; 

) 


C.  DEMODULATION  AND  FILTERING 
1.  DEMODULATE 

This  function  performs  demodulation  of  the  input  carrier  signal  as  well  as  the 
necessary  frequency  shift  when  compensating  for  Doppler.  No  assumption^  are  made 
about  the  sampling  frequency  within  this  program.  It  may  be  set  and  cleared  as  desired. 
The  time  index  n  is  remembered  so  the  repetitive  calls  may  be  made. 


/*  DEMODULATE  performs  complex  demodulation  of  a  sinusoidal  carrier  for 
digital  signal  processing.  Given  an  input  signal,  sampling  frequency, 
and  carrier  frequency  two  outputs  are  produced,  one  ror  the  real  part 
and  one  for  the  imaginary  part.  Tne  user  determines  the  number  of  points 
to  be  processed  in  each  call  to  the  program. 

sin_out  -  Input  signal  and  output  as  imaginary  part  of  the  demodulated 
signal.  Input  data  is  overwritten  and  lost  forever. 

On  initialization  sin_out  is  set  to  the  desired  carrier 
frequency,  (floating  type  pointer.) 

cos_out  -  Output  signal  as  real  part  of  the  demodulated  signal. 

On  Initialization  cos_out  is  set  to  the  desired  sample  frequency, 
(floating  type  pointer.) 

dopp^bin  -  Denotes  the  doppler  frequency  shift  an  the  sampling  frequency 
to  be  used  when  doppler  processing  is  used,  and  is  determined 
externally.  When  no  doppler  processing  is  to  be  used  it  must 
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be  set  to  0. 


n  -  counter  that  maintains  time  increment  for  digit  processing  (static) . 
theta  -  digital  frequency  (static) . 

flag  -  Set  when  system  has  been  initialized  (static) .  Cleared  when 
N  has  been  set  to  zero. 

N  -  the  number  of  data  points  to  be  processed.  Whenever  set  to  zero 
the  system  must  be  reset. 

Reset:  call  with  all  arguments,  N=0. 

Set  : First  call  with  sin_out=fc,  cos_out=fs,  N  is  don't  care. 

Operation  :  Call  with  sin_out  as  input  data,  returns  sin_out  as  imaginary 
part  and  cos_out  as  real  part .  N=  number  of  data  points  to 
process  and  is  the  same  length  as  the  input  data.  */ 

^include  <math.h> 

^define  PI  3.1415926536 

void  demodulate (sin_out, cos_out, dopp_bin, N) 
float  *sin_out,  *cos_out,  dopp_bin; 
int  N; 

{ 

int  i; 

static  int  flag,  n; 
static  double  theta; 

if  (flag  ==  0) 

( 

flag  *>  1; 

theta  =  2  *  PI  * (*sin_out/ (*cos_out) ) ; 

n  =  0; 

return; 

) 

else  if  (  N  ==  0) 

{ 

flag  =  0; 
cheta  =  0.0; 
n  -  0; 
return; 

) 

else 

( 

for(i*=0;  i<N;  i++) 

< 

*cos_out++  •»  *sin_out  *  cos  (n*  (theta  +  dopp_bin) )  ; 

*sin_out++  »  *sin_out  *  sin (n* (theta  +  dopp_bin) ) ; 
n++; 

) 

return; 

) 

) 
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2.  HIGHPASS 


This  program  performs  a  highpass  filtering  of  input  data.  It  is  assumed  to  be 
used  in  conjunction  with  a  lowpass  filter  to  form  a  bandpass  filter.  Processing  in  this 
manner  simplifies  the  processing  by  avoiding  complex  operations.  It  maintains  its  state 
on  subsequent  calls  and  can  be  initialized  with  any  type  filter  up  to  order  20.  It  may 
also  be  cleared  and  reset.  (Note:  All  filtering  programs  are  generic  and  may  be  adapted 
for  any  filtering  application,  lowpass,  highpass,  bandpass  and  bandstop.  The  names 
were  chosen  to  distinguish  them  from  each  other  within  SEQREM.) 


/*  HIGHPASS  is  a  filtering  program  that  operates  using  difference 
equations.  It  is  intended  that  the  filter  be  HR  but  this  is  not 
required,  (i.e.  H(z)  =  B(z)/A(z)  is  a  polynomial  where  the  Bn's  and 
An's  specify  the  coefficients  of  the  difference  equation,  which  must 
of  the  same  length.  If  they  are  not  of  the  same  order  zeros  must  be 
appended. ) 

Initialization:  The  first  call  sets  the  desired  filter  coefficients 
to  be  used  throughout.  No  filtering  is  performed  at  this  stage. 

In  this  case  in_num  specifies  numerator  and  out_den  specifies  the 
denominator  coefficients.  M  is  the  number  of  coefficients,  ordered 
highest  to  lowest.  (NOTE:  out_den(0)  is  assumed  to  be  one  and  is  not 
actually  set  during  the  initialization.) 

Example  1:  H(z)  ■  1/z;  =>  M  =  2,  in_num[0)  =  0,  in_num(l)  =  1, 

out_den[0)  =  1,  out_den[l]  =  0 

Example  2:  H(z)  =  (0.5z*2  +  l)/(zA2  +  0.5z  +  0.6)  =>  M  =  3, 

in_num[0)  =  0.5,  in_num(l)  =  0,  in_num[2)  =  1, 
out_den(0)  =  1,  out_den(l)  =  0.5,  out_den(2)  =0.6 


Filtering:  Subsequent  calls  perform  the  desired  filtering.  In  this 
case  in_num  is  the  input  data,  out_den  is  the  output  data  from  the 
filter,  M  specifies  the  number  of  data  points  being  processed,  and 
may  be  of  any  length  as  long  as  sufficient  space  has  been  allocated 
to  the  input  pointer  and  the  output  pointer. 

Clearing:  To  clear  and  reinitialize  the  filter  after  initializing  once, 
simply  set  M  to  0  and  call  again  as  explained  above. 

Data  Types:  Inputs  are  pointers  to  array's  of  float,  M  is  integer. 

Filter  Order:  Maximum  order  is  25. 

Variables : 

in_num  -  pointer  to  input  numerater  coefficients,  or  to 
input  data  segment. 
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in_den  -  pointer  to  input  denominator  coefficients,  or  to 
output  data  segment.  (NOTE:  input  and  output  segments 
may  be  the  same,  but  input  values  will  be  overwritten.) 

A  -  Denominator  filter  coefficients  (static) . 

B  -  Numerator  filter  coefficients  (static) . 

Y  -  Output  storage  buffer  used  in  recursion  (static) . 

X  -  Input  storage  buffer  used  in  recursion  (static) . 

N  -  Remembers  filter  order  for  computing  output  (static) . 

M  -  Filter  order  during  initialization,  and  number  of  points 
to  be  processed  on  subsequent  calls. 

flag  -  maintains  initialization  status  (static) . 

Other  Filters:  LOWPASSH),  L0WPASS2  () ,  and  L0WPASS3  ()  are  coded 
identically  to  this  program  and  work  in  the  same 
manner.  */ 


{{include  <math.h> 
{{define  MAXORD  25 
{(define  MAX  24 


void  hipass (injnum, out_den,M) 
float  *in_num,  *out_den; 
int  M; 

( 

static  float  A (MAXORD],  B( MAXORD],  Y{ MAXORD),  X[ MAXORD); 
static  int  flag,  N; 
int  i, j; 


/*  Loop  clears  output  und  coefficient  values  for  reuse  in  new  filter.*/ 
if  (M==0) 

( 

for  (i°0;i<MAXORD;i++) 


( 

A  ( i ]  =  0; 
B[i]  =  0; 
Y(i]  =  0; 
X [i]  =  0; 
) 

flag  =0; 
return; 


) 

/*  Loop  sets  coefficients  of  the  desired  filter  on  first  entry  if  flag=0  */ 
/*  Note  that  A { 0 )  is  assumed  1,  since  it  corresponds  to  the  desired  output*/ 
else  if  (flag  ==  0) 


( 

B(0]  •*  in_num(0]; 
for  (i=l;  i<M;i++) 


( 

B(i)  = 
A(i]  ° 
) 

flag  =  1; 
N  ■  M; 


in_num(i)  ; 
out_den (ij ; 
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return; 

) 

/*  Otherwise  perform  filtering  operations.  Output  stored  in  out_den [ ) .  */ 
else 
{ 

for(j=0;  j<M;  j++) 

{ 

X(MAX]  =  in_num(j]; 

/*  Start  filteriing  operation.  */ 

Y [MAX)  =  B[0]  *  X (MAX] ; 
for  (i=l;i<N;i++) 

Y (MAX)  «  Y (MAX)  +  B (i) *X (MAX-i)  -  A (i) *Y(MAX-i) ; 
out_den ( j ]  »  Y (MAX) ; 

/*  Perform  time  shift  of  data  points  stored  in  filter.  */ 
for  (i=0;i<MAX;i++) 

( 

Y(i)  =  Y(i+1); 

X(i)  =  X(i+1 ] ; 

) 

} 

return; 

> 

} 


3.  LOWPASS1 

This  program  performs  a  lowpass  filtering  of  input  data.  It  operates  exactly 
the  same  as  HIGHPASS. 


/*  LOWPASSl  is  a  filtering  program  that  operates  using  difference 
equations.  It  is  intended  that  the  filter  be  HR  but  this  is  not 
required,  (i.e.  H(z)  °  B(z)/A(z)  is  a  polynomial  where  tne  Bn's  and 
An's  specify  the  coefficients  of  the  difference  equation,  which  must 
of  the  same  length.  If  they  are  not  of  the  same  order  zeros  must  be 
appended. ) 

Initialization:  The  first  call  sets  the  desired  filter  coefficients 
to  be  used  throughout.  No  filtering  is  performed  at  this  stage. 

In  this  case  in_num  specifies  numerator  and  out_den  specifies  the 
denominator  coefficients.  M  is  the  number  of  coefficients,  ordered 
highest  to  lowest.  (NOTE:  out_den(0]  is  assumed  to  be  one  and  is  not 
actually  set  during  the  initialization.) 

Example  1:  H(z)  =  1/z;  =>  M  =  2,  in_num(0)  =  0,  in_num(l)  =  1, 

out_den(0)  =  1,  out_den[l)  =  0 

Example  2:  H(z)  =  <0.5zA2  +  l)/(z/>2  +  0.5z  +  0.6)  *■>  M  =  3, 

in_num(0]  =  0.5,  in_num(l)  **  0,  in_num[2)  =  1, 
out_den[0)  =  1,  out_den[l)  *■  0.5,  out_den[2]  =  0.6 


Filtering:  Subsequent  calls  perform  the  desired  filtering.  In  this 
case  in_num  is  the  input  data,  out_den  is  the  output  data  from  the 
filter,  M  specifies  the  number  of  data  points  being  processed,  and 
may  be  of  any  length  as  long  as  sufficient  space  has  been  allocated 
to  the  input  pointer  ana  the  output  pointer. 
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Clearing:  To  clear  and  reinitialize  the  filter  after  initializing  once, 
simply  set  H  to  0  and  call  again  as  explained  above. 

Data  Types:  Inputs  are  pointers  to  array's  of  float,  M  is  integer. 

Filter  Order:  Maximum  order  is  25. 

Variables: 

in_num  -  pointer  to  input  numerater  coefficients,  or  to 
input  data  segment. 

in_den  -  pointer  to  input  denominator  coefficients,  or  to 
output  data  segment.  (NOTE:  input  and  output  segments 
may  be  the  same,  but  input  values  will  be  overwritten.) 

A  -  Denominator  filter  coefficients  (static) . 

B  -  Numerator  filter  coefficients  (static) . 

Y  -  Output  storage  buffer  used  in  recursion  (static) . 

X  -  Input  storage  buffer  used  in  recursion  (static) . 

N  -  Remembers  filter  order  for  computing  output  (static) . 

M  -  Filter  order  during  initialization,  and  number  of  points 
to  be  processed  on  subsequent  calls. 

flag  -  maintains  initialization  status  (static) . 

Other  Filters:  HIGHPASSO,  L0WPASS2 () ,  and  L0WPASS3 ()  are  coded 
identically  to  this  program  and  work  in  the  same 
manner.  */ 

((include  <math.h> 

((define  MAXORD  25 
((define  MAX  24 

void  lowpassl (in_num,out_den,M) 
float  *in_num,  *out_den; 
int  M; 

{ 

static  float  A (MAXORD) ,  B [MAXORD) ,  Y (MAXORD) ,  X(MAXORD); 
static  int  flag,  N; 
int  i, j; 

/*  Loop  clears  output  and  coefficient  values  for  reuse  in  new  filter.*/ 
if  (M==0) 

( 

for  (i“0;i<MAXORD;i+i-) 

( 

A  ( i )  =  0; 

B(i)  =  0; 

Y(i)  =  0; 

X(iJ  =  0; 

) 

flag  =  0; 
return; 

) 

/*  Loop  sets  coefficients  of  the  desired  filter  on  first  entry  if  flag=0  * 
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/*  Note  that  A [ 0 )  is  assumed  1,  since  it  corresponds  to  the  desired  output*/ 
else  if  (flag  ==  0) 

( 

B [ 0  J  =  in_num[0); 
for  (i=l;  i<M;i++) 

( 

B[i)  =  in_num(i); 

A(i)  =  out_den(il; 

) 

flag  =  1; 

N  =  M; 
return; 

} 

/*  Otherwise  perform  filtering  operations.  Output  stored  in  out_den[).  */ 
else 
( 

for  ( j*=0;  j<H;  j++) 

{ 

X  (MAXJ  «=  in_num  [  j ) ; 

/*  Start  filteriing  operation.  */ 

Y (MAX)  =  B [ 0 ]  *  X (MAX) ; 
for  (i=l;i<N;i++) 

Y (MAX)  -  Y (MAXJ  +  B(i] *X (MAX-i)  -  A [i) *Y (MAX-i) ; 
out_den ( j )  =  Y (MAX) ; 

/*  Perform  time  shift  of  data  points  stored  in  filter.  */ 
for  (i=0; i<MAX;i++) 

{ 

Y [i]  «*  Y [ i+1 1 ; 

X[i)  -  X ( i+1 ] ; 

> 

} 

return; 

} 


4.  LOWPASS2 

This  program  performs  a  lowpass  filtering  of  input  data.  It  operates  exactly 
the  same  as  HIGHPASS. 


/*  LOWPASS2  is  a  filtering  program  that  operates  using  difference 
equations.  It  is  intended  that  the  filter  be  IIR  but  this  is  not 
required,  (i.e.  H(z)  =  B(z)/A(z)  is  a  polynomial  where  the  Bn's  and 
An's  specify  the  coefficients  of  the  difference  equation,  which  must 
of  the  same  length.  If  they  are  not  of  the  same  order  zeros  must  be 
appended.) 

Initialization:  The  first  call  sets  the  desired  filter  coefficients 
to  be  used  throughout.  No  filtering  is  performed  at  this  stage. 

In  this  case  in_num  specifies  numerator  and  out_den  specifies  the 
denominator  coefficients.  M  is  the  number  of  coefficients,  ordered 
highest  to  lowest.  (NOTE:  out_den[0)  is  assumed  to  be  one  and  is  not 
actually  set  during  the  initialization.) 

Example  1:  H(z)  =  1/z;  =>  M  =  2,  in_num[0)  »  0,  in_r.um(l)  =  1, 

out_den(0)  «»  1,  out_den(l)  =  0 
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Example  2:  H(z)  =  <0.SzA2  +  l)/(zA2  +  0.5z  +  0.6)  =>  M  =  3, 

in_num{0]  =  0.5,  in_num(l)  =  0,  in_numi2)  =  1, 
out_den(0)  =  1,  out_den(l]  =  0.5,  out_den(2)  =0.6 


Filtering:  Subsequent  calls  perform  the  desired  filtering.  In  this 
case  in_num  is  the  input  data,  out_den  is  the  output  data  from  the 
filter,  M  specifies  the  number  of  data  points  being  processed,  and 
may  be  of  any  length  as  long  as  sufficient  space  has  been  allocated 
to  the  input  pointer  and  the  output  pointer. 

Clearing:  To  clear  and  reinitialize  the  filter  after  initializing  once, 
simply  set  M  to  0  and  call  again  as  explained  above. 

Data  Types:  Inputs  are  pointers  to  array's  of  float,  M  is  integer. 

Filter  Order:  Maximum  order  is  25. 

Variables: 

in_num  -  pointer  to  input  numerator  coefficients,  or  to 
input  data  segment. 

in_den  -  pointer  to  input  denominator  coefficients,  or  to 
output  data  segment.  (NOTE:  input  and  output  segments 
may  be  the  same,  but  input  values  will  be  overwritten.) 

A  -  Denominator  filter  coefficients  (static) . 

B  -  Numerator  filter  coefficients  (static) . 

Y  -  Output  storage  buffer  used  in  recursion  (static) . 

X  -  Input  storage  buffer  used  in  recursion  (static) . 

N  -  Remembers  filter  order  for  computing  output  (static) , 

M  -  ".Iter  order  during  initialization,  and  number  of  points 
to  be  processed  on  subsequent  calls. 

flag  -  maintains  initialization  status  (static) . 

Other  Filters:  HIGHPASSO,  LOWPASSK),  and  LOWPASS3  ()  are  coded 
identically  to  this  program  and  work  in  the  same 
manner.  */ 

((include  <math.h> 

((define  MAXORD  25 
((define  MAX  24 

void  lowpass2 (in_num,out_den,M) 
float  *in_num,  *out_den; 
int  M; 

{ 

static  float  A (MAXORD),  B (MAXORD),  Y (MAXORD),  X(MAXORD); 
static  int  flag,  N; 
int  i, j; 

/*  Loop  clears  output  and  coefficient  values  for  reuse  in  new  filter.*/ 
if  (M«=0) 

( 
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for  (i=0;i<MAXORD;i++) 

{ 

A(i]  =  0; 

B[i)  -  0; 

Y{i]  -  0; 

X[i)  «=  0; 

) 

flag  ■  0; 
return; 

} 

/*  Loop  sets  coefficients  of  the  desired  filter  on  first  entry  if  flag=0  */ 
/*  Note  that  A ( 0 J  is  assumed  1,  since  it  corresponds  to  the  desired  output*/ 
else  if  (flag  ==  0) 

{ 

B[0)  =  in_num[0); 
for  (i=l;  i<M;i++) 

( 

B[i]  =  in_num(i); 

A[i]  =  out_den(i); 

) 

flag  ■  1; 

N  ■=*  M; 
return; 

} 

/*  Otherwise  perform  filtering  operations.  Output  stored  in  out_den[].  */ 
else 
( 

for(j=0;  j<M;  j++) 

( 

X{MAX]  =  in_num(j); 

/*  Start  filteriing  operation.  */ 

Y(MAX)  =  B { 0 }  *  X (MAX) ; 
for  (i»l;i<N;i++) 

Y(MAX)  =  Y (MAXJ  +  B(i] *X (MAX-i]  -  AJi) *Y (MAX-i) ; 
out_den ( j ]  =  Y [MAX] ; 

/*  Perform  time  shift  of  data  points  stored  in  filter.  */ 
for  (i“0;i<MAX;i++) 

( 

Y [ i J  =  Y [i+1 ) ; 

X(i]  ■=  X C i+1 )  ; 

) 

) 

return; 

) 

} 


5.  LOWPASS3 

This  program  performs  a  lowpass  filtering  of  input  data.  It  operates  exactly 
the  same  as  HIGHPASS. 


/*  L0WPASS3  is  a  filtering  program  that  operates  using  difference 
equations.  It  is  intended  that  the  filter  be  HR  but  this  is  not 
required,  (i.e.  H(z)  =  B(z)/A(z)  is  a  polynomial  where  the  Bn's  and 
An's  specify  the  coefficients  of  the  difference  equation,  which  must 
of  the  same  length.  If  they  are  not  of  the  same  order  zeros  must  be 
appended . ) 
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Initialization:  The  first  call  sets  the  desired  filter  coefficients 
to  be  used  throughout.  No  filtering  is  performed  at  this  stage. 

In  this  case  in_num  specifies  numerator  and  out_den  specifies  the 
denominator  coefficients.  M  is  the  number  of  coefficients,  ordered 
highest  to  lowest.  (NOTE:  out_den(0)  is  assumed  to  be  one  and  is  not 
actually  set  during  the  initialization.) 

Example  1:  H(z)  =  1/z;  =>  M  ■=  2,  in_num[0)  =  0,  in_num[l)  ■=  1, 

out_den(0]  =  1,  out_den(l)  =  0 

Example  2:  H(z)  =  <0.5z*2  +  l)/(zA2  +  0.5z  +  0.6)  =>  M  =  3, 

in_num(0]  =  0.5,  in_num(l]  =  0,  in_num(2)  ■=  1, 
out_den[0)  =  1,  out_den(l)  =  0.5,  out_den(2)  =  0.6 


Filtering:  Subsequent  calls  perform  the  desired  filtering.  In  this 
case  in_num  is  the  input  data,  out_den  is  the  output  data  from  the 
filter,  M  specifies  the  number  of  data  points  being  processed,  and 
may  be  of  any  length  as  long  as  sufficient  space  has  been  allocated 
to  the  input  pointer  and  the  output  pointer. 

Clearing:  To  clear  and  reinitialize  the  filter  after  initializing  once, 
simply  set  M  to  0  and  call  again  as  explained  above. 

Data  Types:  Inputs  are  pointers  to  array's  of  float,  M  is  integer. 

Filter  Order:  Maximum  order  is  25. 

Variables: 

in_num  -  pointer  to  input  numerater  coefficients,  or  to 
input  data  segment. 

in_den  -  pointer  to  input  denominator  coefficients,  or  to 
output  data  segment.  (NOTE:  input  and  output  segments 
may  be  the  same,  but  input  values  will  be  overwritten.) 

A  -  Denominator  filter  coefficients  (static) . 

B  -  Numerator  filter  coefficients  (static) . 

Y  -  Output  storage  buffer  used  in  recursion  (static) . 

X  -  Input  storage  buffer  used  in  recursion  (static) . 

N  -  Remembers  filter  order  for  computing  output  (static) . 

M  -  Filter  order  during  initialization,  and  number  of  points 
to  be  processed  on  subsequent  calls. 

flag  -  maintains  initialization  status  (static) . 

Other  Filters:  HIGHPASSO,  LOWPASSK),  and  LOWPASS2  ( )  are  coded 
identically  to  this  program  and  work  in  the  same 
manner.  */ 

^include  <math.h> 
fdefine  MAXORD  25 
idefine  MAX  24 
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void  lowpass3 (in_num,out_den,M) 
float  *in_num,  *out_den; 
int  M; 

{ 

static  float  A(MAXORD),  B(MAXORD) ,  Y (MAXORD) ,  X(MAXORD); 
static  int  flag,  N; 
int  i, j; 

/*  Loop  clears  output  and  coefficient  values  for  reuse  in  new  filter.*/ 
if  (M==0) 

( 

for  (i=0;i<MAXORD;i++) 

{ 

A[i]  =  0; 

B[i)  =  0; 

Y[iJ  =  0; 

X  (ij  =  0; 

} 

flag  =  0; 
return? 

} 

/*  Loop  sets  coefficients  of  the  desired  filter  on  first  entry  if  flag=0  */ 
/*  Note  that  A [ 0 ]  is  assumed  1,  since  it  corresponds  to  the  desired  output*/ 
else  if  (flag  ==  0) 

{ 

B [ 0 )  =  in_num[0); 
for  (i=l;  i<M;i++) 

{ 

B{i)  =  in_num[i]; 

Ali)  =  out  denlij; 

} 

flag  “  1; 

N  =  M; 
return; 

} 

/*  Otherwise  perform  filtering  operations.  Output  stored  in  out_den().  */ 
else 
{ 

for ( j=0;  j<M;  j++) 

( 

X(MAX)  =  in_num(jl; 

/*  Start  filteriing  operation.  */ 

Y [MAX]  =  B ( 0 ]  *  X(MAX); 
for  (i=l;i<N;i++) 

Y [MAX)  =  Y(MAX)  +  B(i] *X(KAX-i]  -  A(i] *Y[HAX-i) ; 
out_den ( j  J  =  Y (MAX) ; 

/*  Perform  time  shift  of  data  points  stored  in  filter.  */ 
for  (i=0;i<MAX;i++) 

< 

Y ( i J  =  Y (i+1)  ; 

X ( i )  =  X(i+1) ; 

) 

} 

return; 

) 
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D.  FAST  HADAMARD  TRANSFORM  (FHT) 

1.  HADAMARD 

This  function  performs  the  FHT  of  an  input  data  vector  in  place.  A  power  of 
two  length  is  required.  Computation  is  identical  to  the  FFT  butterfly  process  with  the 
complex  mul  jplies  set  to  one. 


/*  HADAMARD  performs  the  Fast  Hadamard  Transform  (FHT)  on  the  input  data. 
The  transform  is  performed  in  place  and  assumes  no  interleaving. 
Interleaving  must  be  performed  externally. 

The  FHT  performs  the  same  operations  as  the  FFT  with  the  complex 
multiplies  set  to  +  or  -  1 .  The  input  and  output  data  vector 
must  be  first  ordered  and  then  reordered  after  exit. 

Variables: 

degree  -  Provides  the  degree  of  the  polynomial  being  worked  with 
which  gives  the  length  of  the  M-Sequence. 

din  -  Pointer  to  the  input  data  to  be  transformed.  On  exit  it  holds 
the  transformed  data. 

m,n  -  counters. 

itop, ispace, iwidth, ibot  -  common  place  holders  used  performing 
the  butterfly  operations  of  the  FFT.  */ 


einclude  <math.h> 

void  hadamard (degree, din) 
float  *din; 
int  degree; 

( 

int  m, n, itop, ispace, iwidth, ibot; 
float  temp; 

unsigned  int  one,  size; 
one  “1; 

size  =  (one«degree)  ; 
for  (  m=l;  m<“degree;m++) 

( 

ispace  =  (one«m)  ; 
iwidth  =  (one«(m-l)); 
for  (n»0;n<iwidth;n++) 

for(  itop**n;itop<*>  (size-2)  ;itop+=ispace) 

( 

ibot  =  itop  +  iwidth; 
temp  »  din (ibot); 
din(ibot)  =  din (itop)  -  temp; 
din  (itop)  •>  din(itop)  +  temp; 

) 
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} 

} 

E.  MAGNITUDE  AND  PHASE 
1.  MAG  PHASE 

This  function  simply  computes  the  magnitude  and  phase  of  the  demodulates 
after  either  coherent  averaging  or  directly  after  the  FBT  and  the  appropriate 
permutations. 


/*  MAG_PHASE  computes  the  magnitude  and  phase  of  the  input  data,  and 
returns  the  results  in  the  pointers  mag  and  phase. 

realpart  -  Pointer  for  the  real  part  of  the  input  data  of  length  N. 

impart  -  Pointer  for  the  imaginary  part  of  the  input  data  of  length  N. 

mag  -  Pointer  to  data  segment  to  hold  the  magnitude  of  the  data  and 
is  also  of  length  N. 

phase  -  Pointer  for  data  segment  to  hold  the  phase  of  the  data  and  is 
also  of  length  N.  (Degrees) 

N  -  Integer  indicating  the  size  of  the  input  data  segment.  */ 

linclude  <math.h> 

Sdefine  PI  3.1-515927 

void  mag_phase (realpart, impart, mag, phase, N) 
float  *realpart,  *impart,  *mag,  *phase; 
int  N; 

< 

int  i; 

for  (i*»0;i<N;i++) 

( 

mag[i)  =  sqrt (realpart (ij  *  realpart [ij  +  impart (i)  *  impart (il); 
if  (realpart (i1  ==  0.0) 

{ 

phased)  =  0.0; 

) 

else 

phased)  -  atan  (impart  (i)  /  realpart(ij)  *  (180.0/PI); 

) 

) 
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F.  UTILITY  PROGRAMS 
1.  SET  FILTER 

This  is  a  handy  MATLAB  program  for  finding  the  coefficients  for  a  tenth 
order  Butterworth  bandpass  filter  and  a  fifth  order  Chebychev  lowpass  filter,  which 
stores  them  in  a  file  formatted  for  use  with  GET_FLT_CGEF.  Sampling  frequency 
and  cutoff  frequencies  are  prompted  from  the  user. 


%  function  set_filter  Set_filter  computes  the  filter  coefficients  for 
%  the  bandpass  filter  (high  and  lowpass  filters) 

%  and  for  the  lowpass  filter  used  in  the  sequence 

%  removal  program.  The  bandpass  region  uses  a 

%  10th  order  BUTTERWORTH  l  Iter  and  the  lowpass  region 

%  uses  a  5th  order  CHEBYCHEV  filter  with  the  results 

%  stored  in  the  user  specified  file  in  the  format 

%  required  by  the  GET_FLT_CQSF  function.  The  user  is 

%  queried  for  cutoff  frequencies  and  sampling 

%  frequency. 

*  function  set_fiiter 

fprintf ( '\n  WARNING!:;!!!  \n’>; 

fprintf ('\n  File  to  hold  results  should  be  non-existarit  prior  to  run»ing\n') ; 
fprintf ('  this  routine  because  the  results  are  appended  to  any  pre-existing\n ') ; 
fprintf ('  data.Xn'); 

filename  =  input ('Enter  filename  to  store  results:  ’.'s'); 

fs  *»  inoutCEr.ter  the  Sampling  frequency  to  be  used:  ■); 

fc  «*  input(  'Enter  the  Cut  Off  Frequency  for  the  Low  Pass  Filter:  '); 

fprintf('\n  Enter  the  Cut  Off  Frequencies  for  the  BP  filter. \n'); 

fl  =  input (' Lower :  '); 

f h  =  input ( ' Upper :  ■ ) ; 

fs  =  fs/2; 

if  (fl  >  fs)  I  (fc  >  fs)  I  (fh  >  fs) 

error ('Invalid  Frequencies.  Cutoff  Frequencies  Must  be  <=  Fs/2'); 

end 

%  Finds  the  Coefficients  for  the  High  Pass  Filter  (BP  low  end) 
wn  *  fl/fs; 

[bl.al]  «  butrer (10,wn, 'high') ; 

%  Finds  the  Coefficients  for  the  Low  Pass  Filter  (BP  hi  end) 
wn  »  fh/fs; 

(b2,a2)  -  butter (10, wn) ; 

%  Finds  the  Coefficients  for  the  Low  Pass  Filter. 

•  wn  “  fc/fs; 

(b3,a3j  -  chebyl(5,0.5,wn) ; 
for  1*1:11 

fprintf (filename, '%e  %e\n',bl (i) ,al (i) ) ; 

‘  end 

for  1*1:11 
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fprintf (filename, '%e  %e\n’ ,  b2 (i) , a2 (i) ) ; 

end 

for  1*1:6 

fprintf (filename, 1 %e  %e\n ' ,b3 (i) ,a3 (i) ) ; 

end 

2.  APLOT 

This  program  was  used  to  generate  the  plots  in  Appendix  A  with  NGAR 
graphics.  It  is  data  specific  because  of  the  nature  of  NCAR  graphics  when  used  for 
surface  plotting  routines. 


C  This  is  a  program  to  plot  the  results  of  the  program  SEQ_REM 
C  in  successive  plots  for  each  doppler  searched. 

INTEGER  X,Y 

PARAMETER  (X=510,  Y=30) 

CHARACTER*32  IFILE, TEMP, SEQLEN 

CHARACTER*20  DOPPLER 

INTEGER  I, J, LENGTH, ISZ 

REAL  MAG (X, Y) , PH (X, Y) , XPOS, YPOS, ZPOS 

REAL  WORK (2*X*Y+X+Y) 

REAL  TS, ANGH, ANGV, ZMAX 

COMMON  /SRFIP1/  IFR, ISTP, IROTS, IDRX, IDRY, IDRZ, IOPPER,  ISKIRT, 
1  NCLA, THETA, HSKIRT, CHI , CLO, CINC, ISPVAL 

35  FORMAT (A) 

40  FORMAT  (A20) 

IFR  =  0 
IDRZ  =  1 

WRITE**,*) 'THIS  PLOTS  THE  RESULTS  FROM  SEQREM  PROGRAM' 

WRITE  <*,*)-' ENTER  THE  INPUT  FILE:  ' 

READ  (*,40)  IFILE 

OPEN (33, FILE=IFILE, STATUS* 'OLD ' ) 

ANGH  =  -60. 

ANGV  =  30. 

ZMAX  =  0. 

CALL  GOPKS (6, ISZ) 

CALL  GOPWK( 1,2,1) 

CALL  GACWK(l) 

READ (33, 35) TEMP 

50  READ (UNIT=33,FMT=*) DOPPLER 

READ (UNIT=33, FMT=*) SEQLEN 
READ (UNIT=33, FMT=*) LENGTH, TS 
J  =  1 

60  J  -  J+l 

DO  80  1=1, LENGTH 

READ <UNIT=33, FMT=*, ERR=100, END=200) MAG (I,  J) , PH (I,  J) 

IF  (MAG (I, J)  .GT.  ZMAX)  THEN 
ZMAX  =  MAG (I, J) 

ENDIF 

80  CONTINUE 
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GO  TO  60 


100  J  =  J-l 

CALL  GSELNT  (0) 

CALL  GSTXAL(2, 3) 

CALL  GSCHH ( . 02) 

CALL  GTX(0.5, 0.9, 'ARRIVAL  TIME  PLOT  ') 

CALL  GSTXAL (2, 3) 

CALL  GTX (0.5, 0.85, DOPPLER) 

CALL  EZSRFC (MAG, LENGTH, J, ANGH, ANGV, WORK) 

XPOS  =  0. 

YPOS  =  -0.25*ZMAX 
ZPOS  =  0.0 

CALL  PWRZS (XPOS, YPOS, ZPOS, SEQLEN, 31, 25,1,3,0) 
CALL  FRAME 
ZMAX  “  0. 

GO  TO  50 

200  J  =  J-l 

CALL  GSELNT (0) 

CALL  GSTXAL (2, 3) 

CALL  GSCHH (.02) 

CALL  GTX(0. 5, 0.9, 'ARRIVAL  TIME  PLOT  ') 

CALL  GSTXAL (2, 3) 

CALL  GTX (0.5, 0.85, DOPPLER) 

CALL  EZSRFC (MAG, LENGTH, J, ANGH, ANGV, WORK) 

XPOS  •*  0. 

YPOS  *=  -0.25*ZMAX 
ZPOS  =  0.0 

CALL  PWRZS (XPOS, YPOS, ZPOS, SEQLEN, 31, 25, 1, 3, 0) 
CALL  FRAME 
CLOSE (33) 

CALL  GDAWK(l) 

CALL  GCLWK(l) 

CALL  GCLKS 

WRITE (* , *) 'FINISHED' 

STOP 

END 


89 


3.  DOPPLOT 

This  program  was  used  to  generate  the  arrival  time  plot,  Figure  3.1,  with 
NCAR  graphics.  It  is  data  specific  because  of  the  nature  of  NCAR  graphics  when  used 
for  surface  plotting  routines. 


C  This  is  a  program  to  plot  the  results  of  the  program  SEQ_REM 
C  making  a  doppler  versus  time  plot  for  each  doppler  search. 

INTEGER  X, Y 

PARAMETER  (X=255,  Y=23) 

CHARACTER*32  IFILE, TEMP, SEQLEN 
CHARACTER*20  DOPPLER 
INTEGER  I, J, LENGTH, ISZ 

REAL  MAG  (X,  Y)  ,  PH  (X,  Y)  ,  XPOS,  YPOS,  ZPOS,  TEMPI  (X,  Y)  ,  TEMP2  (X,  Y) 

REAL  W0RK(2*X*Y+X+Y) 

REAL  TS, ANGH, ANGV, ZMAX 

COMMON  /SRFIP1/  IFR, ISTP, IROTS, IDRX, IDRY, IDRZ, IOPPER,  ISKIRT, 

1  NCLA, THETA, HSKIRT, CHI, CLO, CINC, ISPVAL 

35  FORMAT (A) 

40  FORMAT  (A20) 

IFR  =  0 
IDRZ  «  1 

WRITE(*,*) ‘THIS  PLOTS  THE  RESULTS  FROM  SEQREM  PROGRAM  IN  DOPPLER' 
WRITEt*,*)  'VS. TIME  FORMAT' 

WRITE (*,*)  'ENTER  THE  INPUT  FILE:  ’ 

READ  (*,40)  IFILE 

OPEN (33, FILE=IFILE, STATUS= 'OLD ' ) 

ANGH  =  -60. 

ANGV  =  30. 

ZMAX  «*  0. 

CALL  GOPKS(6, ISZ) 

CALL  GOPWK( 1,2,1) 

CALL  GACWK(l) 

READ (33, 35) TEMP 
J  -  1 

50  READ (UNIT=33,FMT=*) DOPPLER 
READ (UNIT=33, FMT=») SEQLEN 
READ(UNIT=33, FMT=*) LENGTH,  TS 
K  =  1 

DO  53  1=1, LENGTH 
READ  ( 33,  * )  TEMPI  (I,  K) ,  TEMP2  ( I ,  K) 

53  CONTINUE 

60  J  =  J+l 

DO  80  1=1, LENGTH 

READ (UNIT=33, FMT=*, ERR=50, END=200) MAG (I, J) , PH (I, J) 

IF  (MAG(I, J)  .GT.  ZMAX)  THEN 
ZMAX  =  MAG (I, J) 

ENDIF 

80  CONTINUE 

82  CONTINUE 
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DO  85  1=1, LENGTH 

READ (UNIT=33, FMT=*, ERR=50, END=200) TEMPI (I, K) , TEMP2 (I, K) 
85  CONTINUE 
GO  TO  82 


200  J  =  J-l 

CALL  GSELNT (0) 

CALL  GSTXAL(2, 3) 

CALL  GSCHH { . 02) 

CALL  GTX(0.5, 0.9, 'ARRIVAL  TIME  VS  DOPPLER') 

CALL  GSTXAL(2,  3) 

CALL  GTX (0.5, 0.85, '-5  TO  5  KNOTS  DOPPLER') 

CALL  EZSRFC (MAG, LENGTH, J, ANGH, ANGV, WORK) 

XPOS  =  0. 

YPOS  =  -0.25*ZMAX 
ZPOS  =0.0 

CALL  PWRZS (XPOS, YPOS, ZPOS, 'ARRIVAL  TIME  22.36  SEC. ' , 31, 25, 1, 3, 0) 
CALL  FRAME 
CLOSE (33) 

CALL  GDAWK(l) 

CALL  GCLWK(l) 

CALL  GCLKS 

WRITE (*, *) 'FINISHED' 

STOP 

END 


4.  MGEN 

MGEN  generates  a  file  with  a  m-sequence  in  the  forward  direction  in  column 
form,  given  a  primitive  polynomial.  The  sequence  is  in  (1,-1)  form.  A  second  file  is 
also  created,  which  holds  one  period  of  the  register  states.  This  can  be  useful  for 
debugging  and  as  an  aid  for  understanding  m-sequences. 


/*  This  is  a  generator  routine  to  generate  the  output  from  the  least 
significant  register  in  a  M-Sequence  shifter  register  generator. 

Output  is  in  the  form  of  -1  and  1  where  the  mapping  is  (-1,1)  -> 

(1,0).  Register  states  are  also  stored  to  a  specified  file  in  decimal 
form. 

Input:  law  =  the  polynomial  law  specifying  the  M-Sequence  to  be 
generated. 

initial  =  the  initial  load  placed  in  the  Shift  Registers. 

Output:  output  is  a  column  format  of  -1  and  i's  in  floating  format 

form  for  ease  of  use.  output  2  is  a  column  of  numbers  in  decimal 
representing  the  register  states  over  one  period. 

Initiation:  use  'MGEN' 

*/ 

^include  <malloc.h> 

Sinclude  <stdio.h> 

I include  <math.h> 
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main  () 

{ 

float  *mout; 

unsigned  int  temp,  initial,  law,  seq_len; 
int  i; 

char  filename (20) ; 

FILE  *fpl, *fp2; 

printf ("\nEnter  the  Polynomial  Law  in  Octal  Form:  ") ; 
scanf ("%o", slaw) ; 

printf ("\nEnter  the  Initial  Register  Load  in  Octal (Not  2ero!):  "); 
scanf ("%o", sinitial) ; 
if  (initial==0) 

{ 

printf ("\n0  is  Invalid  Register  Load.  Aborting! \n") ; 
exit  (1) ; 

) 

temp  =  law; 
seq_len  =  1; 
while  (temp»=l) 
seq_len<<=l; 
seq_len — ; 

if  ( (mout  =  (float  *)malloc(seq_len*sizeof (mout) ) )==N0LL) 

( 

printf (“\nCannot  Allocate  Memory  to  Output  VariableXn") ; 
exit  (1) ; 

) 

printf ("NnEnter  file  to  store  M-sequence.  (20  Characters  max."); 
scanf ("%s", filename) ; 
fpl=  fopen (filename, "w") ; 

printf ("\nEnter  file  to  store  register  states.  (20  Characters  max.") 
scanf ("%s", filename) ; 
fp2“  fopen(filename,"w"); 
for  (i=0;i<sec;>_len;i++) 

( 

fprintf (fp2, "%u\n", initial) ; 
if  (initialSl) 

{ 

mout[iJ  =  -1.0; 

initial  =  (initial^law) »1; 

) 

else 

( 

mout ( i )  =  1.0; 
initial»=l; 

) 

) 

f close (fp2) ; 
f or (i=0; i<seq_len; i++) 

fprintf (fpl, "%f\n", *mout++) ; 
fclose(fpl) ; 

) 
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5.  MAKEFILE 


This  is  a  file  which  makes  compiling  and  linking  of  C  programs  easier  on 
UNIX  based  systems.  Files  in  the  list  are  compiled,  if  not  already,  and  are  linked 
automatically  with  the  MAKE  command.  The  executable  code  is  stored  in  the  file 
following  the  -o  flag,  and  -lm  includes  the  necessary  library  routines. 

^makefile  for  thesis  work. 

OBJECTS  =  seqrem.o  rev_had.o  init_had.o  fwd_had.o  get_flt_coef .o  low_filterl.o  low_filter2 .o 
low_filter3.o  hi_filter.o  demod.o  hadamard.o  mag_phase.o 

SOURCES  =  seqrem.c  set_had.c  init_had.c  fwd_had.c  get_flt_coef  .c  low_filterl. c  lov>_j:ilter2  .c 
low_filter3.c  hi_f liter. c  demod.c  hadamard.o  mag_phase.c 

LIBS  =  -lm 

CFLAGS  =  -O 

seqrem:  $ (OBJECTS) 

cc  $  (CFLAGS)  S (OBJECTS)  -o  seqrem  $(LIBS) 

$ (OBJECTS):  macrofile. h 


6.  EXAMPLE  FILTER  COEFFICIENT  FILE 

This  shows  the  structure  of  the  format  used  to  store  the  filter  coefficients  used 
in  the  various  filtering  programs.  The  numerator  coefficients  constitute  the  first 
column  and  the  denominator  the  second.  The  first  eleven  rows  are  for  the  first  tenth 
order  filter,  the  next  the  second,  and  finally  the  last  six  make  up  the  fifth  order  lowpass 
filter. 


93 


6.127443e-03 

-6.127443e-02 

2.757350e-01 

-7.352932e-01 

1.286763e+00 

-1.544116e+00 

1.286763e+00 

-7.352932e-01 

2 . 757350e-01 

-6.127443e-02 

6.127443e-03 

6. 127443e-03 

6. 127443e-02 

2.757350e-01 

7 . 352932e-'01 

1.286763e+00 

1.544116e+00 

1.286763e+00 

7.352932e-01 

2 .  757350e-01 

6. 127443e-02 

6.127443e-03 

1.438164e-05 

7.190822e-05 

1 . 438164e-04 

1 . 438164e-04 

7 . 190822e-05 

1 . 438164e-0£ 


1.000000e+00 
-9. 960099e-01 
1 . 759559e*00 
-1 . 112098e+00 
8 . 747430e-01 
-3 . 474470e-01 
1 . 441634e-01 
-3.307116e-02 
6. 691890e-03 
-6.798897e-04 
3.9147X0e-05 
1 . 000000e+00 
9. 960099e-01 
1.759559e+00 
1.112098e+00 
8.747430e-01 
3.474470e-01 
1.441634e-01 
3.307116e-02 
6. 691890e-03 
6.798897e-04 
3 . 9x4710e-05 
1.000000e+00 
-4.515326e+00 
8.278174e+00 
-7 . 695624e+00 
3 . 625292e+00 
-6. 920557e-01 
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